key: cord-0806278-5igf77tk
authors: Zhang, Wenjing; Jang, Sophia; Jonsson, Colleen B; Allen, Linda J S
title: Models of cytokine dynamics in the inflammatory response of viral zoonotic infectious diseases
date: 2018-06-29
journal: Math Med Biol
DOI: 10.1093/imammb/dqy009
sha: 99d431d905c498c0627d077c5c3a69bdb729048b
doc_id: 806278
cord_uid: 5igf77tk

Inflammatory responses to an infection from a zoonotic pathogen, such as avian influenza viruses, hantaviruses and some coronaviruses, are distinctly different in their natural reservoir versus human host. While not as well studied in the natural reservoirs, the pro-inflammatory response and viral replication appear controlled and show no obvious pathology. In contrast, infection in humans results in an initial high viral load marked by an aggressive pro-inflammatory response known as a cytokine storm. The key difference in the course of the infection between the reservoir and human host is the inflammatory response. In this investigation, we apply a simple two-component differential equation model for pro-inflammatory and anti-inflammatory responses and a detailed mathematical analysis to identify specific regions in parameter space for single stable endemic equilibrium, bistability or periodic solutions. The extensions of the deterministic model to two stochastic models account for variability in responses seen at the cell (local) or tissue (global) levels. Numerical solutions of the stochastic models exhibit outcomes that are typical of a chronic infection in the natural reservoir or a cytokine storm in human infection. In the chronic infection, occasional flare-ups between high and low responses occur when model parameters are in a region of bistability or periodic solutions. The cytokine storm with a vigorous pro-inflammatory response and less vigorous anti-inflammatory response occurs in the parameter region for a single stable endemic equilibrium with a strong pro-inflammatory response. The results of the model analyses and the simulations are interpreted in terms of the functional role of the cytokines and the inflammatory responses seen in infection of the natural reservoir or of the human host.

Many human infectious diseases originate as zoonoses through the cross-species transmission from wildlife to humans. Often they are ribonucleic acid (RNA) viral zoonoses, e.g. highly pathogenic avian influenza viruses (HPAIV), hantaviruses (HTV) and severe acute respiratory syndrome coronavirus (SARS-CoV) (Jones et al., 2008; Lloyd-Smith et al., 2009) . They pose a significant threat to public health worldwide and are classified as emerging infectious diseases (Jones et al., 2008; Wang & 270 W. ZHANG ET AL. Crameri, 2014) . The spillover of these viruses to humans m ay result in severe, life-threatening disease, e.g. highly pathogenic avian influenza (HPAI), hantavirus cardiopulmonary syndrome (HCPS), hemorrhagic fever with renal syndrome (HFRS) and severe acute respiratory syndrome (SARS).

Despite major differences among HPAIV, HTV and SARS-CoV, in their coding and replication strategies and the activation of specific pro-and anti-inflammatory cytokines, similar inflammatory responses are seen from infection with these viruses. The site of viral infection and disease associated with HPAI, HCPS, HFRS and SARS is the lower respiratory tract. It is generally assumed that the cytokine response drives the pathology rather than the viral load (Liu et al., 2016) . The major target cells for HPAIV and SARS-CoV replication are type II aveolar epithelial cells (Weinheimer et al., 2012; Qian et al., 2013) and for HTV they are microvascular endothelial (Mackow et al., 2013) . Other cells, such as macrophages and dendritic cells, may become infected and/or release pro-inflammatory cytokines that in turn affect their own cells (autocrine signalling) or other cells (paracrine signalling) to release additional cytokines (He et al., 2006) . The clinical syndrome is characterized by acute lung injury, respiratory failure and a cytokine storm (Chan et al., 2005; Tisoncik et al., 2012) . The cytokineinduced pro-inflammatory response plays a significant role in the pathology of these diseases (Frieman et al., 2008; Peiris et al., 2009 ) and leads to a system-wide cascade effect known as a cytokine storm (Tisoncik et al., 2012) .

In viral zoonoses, the term 'reservoir' implies that the virus is maintained in nature by a particular host. The natural reservoir for HTV is wild rodents, for HPAIV, it is wild aquatic birds (World Health Organization, 2005; Wang et al., 2016) and for SARS-CoV, it is bats (Li et al., 2005; Hu et al., 2015) . Humans are not the natural reservoir and therefore, they only become infected through direct or indirect contact (secondary host or environmental source) with the natural reservoir. Transmission occurs through direct contact with the reservoir or indirectly through the environment such as inhalation of virus shed through excreta which may be from the natural reservoir or an intermediate carrier (domestic poultry in avian influenza or palm civet in SARS-CoV). For SARS-CoV, the virus also spreads via direct person-to-person contact.

The virus often persists for the life of the animal in the natural reservoir (Meyer and Schmaljohn, 2000; World Health Organization, 2005; Vorou, 2016) . In the case of Andes virus (ANDV), according to Jonsson et al. (2010) at least 21 HTV have been associated with human illness. Each of these HTV has a unique rodent reservoir, e.g. Sin Nombre virus's reservoir is the deer mouse, ANDV's reservoir is the long-tailed pygmy rice rat and Black Creek Canal virus's reservoir is the cotton rat. In avian influenza viruses, wild waterfowl, gulls and shorebirds are the natural reservoirs (World Health Organization, 2005) . Horseshoe bats in China have been identified as reservoirs for a novel coronavirus related to SARS-CoV (Hu et al., 2015) .

There are distinct differences in the immune response to infection in the natural reservoir versus spillover to humans. In the natural reservoir for HTV, there is little to no pro-inflammatory response during the early stages in the lung; but in later stages of the infection, anti-inflammatory cytokine levels and tumour growth factor (TGF-β) increase (Schountz et al., 2012) . These cytokines can influence the maintenance or differentiation of other cells, e.g. regulatory T cells (Tregs). Tregs may act indirectly allowing a persistent infection by limiting T helper cell responses or by modulating antigen-presenting cell function or may act directly by cell-to-cell contact (Easterbrook et al., 2007; Easterbrook and Klein, 2008) . HPAI H5N1 infection in ducks, the natural reservoir, versus chickens shows different patterns of expression of the interferon-induced transmembrane protein (IFITM) gene family; IFITM1, 2 and 3 genes are upregulated in response to HPAIV in ducks but not in chickens (Smith et al., 2015) . Bats, the natural reservoir for SARS-CoV, show no clinical signs of disease, and the mechanisms that maintain the virus are not known (Wynne & Wang, 2013) .

In contrast to the natural reservoir, the innate immune response in human infection is distinctly different. It has been shown that human infection with the viruses associated with HCPS, HPAI or SARS results in upregulation of the pro-inflammatory cytokines, a major contributing factor to increased disease severity (He et al., 2006; Szretter et al., 2007; Frieman et al., 2008; Peiris et al., 2009; Safronetz et al., 2011) . In HPAI H5N1 virus infection of human alveolar epithelial cells, IL-6, IL-8 and IFN-β were upregulated, while infection of primary human macrophages induced TNF-α, IL-1, IFN-α, IFN-β and IL-1β (Szretter et al., 2007) . In SARS-CoV infected cells, high levels of MCP-1 and TGF-β1 and intermediate levels of TNF-α, IL-1β and IL-6 were detected (He et al., 2006) . Soluble cytokine receptors such as sTNFR1 or SIL-1RII may also play a role in the regulation of the proinflammatory response by competing with TNF or IL1 and prevent their binding of their receptor. Herein, we define this as an anti-inflammatory activity, (Heaney & Golde, 1996) and we focus on this effect in our modelling.

Key differences of the innate immune response, namely, the magnitude of the pro-and antiinflammatory responses, in natural reservoirs versus human hosts have not been modelled. The major goal of this investigation is to propose a simple modelling framework with pro-inflammatory and anti-inflammatory cytokines to account for the differences between the natural reservoir, where the interaction between pro-inflammatory and anti-inflammatory cytokines maintains a low level of cytokines (characteristic of a chronic but mild persistent infection) and the human host, where a vigorous pro-inflammatory response leads to a cytokine storm. In addition, various models illustrate that alternative pathways of the innate immune response may lead to similar outcomes. We apply a model formulated by Baker et al. (2012) that includes two variables, a system of ordinary differential equations (ODEs), one for the concentration of pro-inflammatory cytokines and the second one for anti-inflammatory cytokines. Hereafter, this model is referred to as the Baker model. This twocompartment model is a non-mechanistic approximation to the real system. This simple model with just a few parameters enables us to understand some basic interactions between the cytokines. The mechanisms responsible for the cytokine storm are much more complex. With time series data from carefully designed laboratory experiments, specific pathways and genes/proteins associated with proinflammatory and anti-inflammatory responses can be identified and used to test some of the outcomes of this simple model. The interaction of these two broad categories of cytokines is posited to play important roles in inflammatory immune responses (Seymour & Henderson, 2001) for both acute infectious diseases (Henderson et al., 1998) and chronic autoimmune diseases (McInnes & Schett, 2007) .

The Baker et al. (2012) study focuses on a chronic inflammatory disease, rheumatoid arthritis. Instead, we apply their model to describe two different responses, one that mimics the controlled response in the natural reservoir and the second one that mimics the out-of-control, pro-inflammatory response in human infection. We focus primarily on the relation between the parameter that upregulates the pro-inflammatory response and three other parameters: representing clearance ratio, upregulation of anti-inflammatory response and background production of pro-inflammatory cytokines. We concentrate on their functional role rather than identifying specific pro-inflammatory or anti-inflammatory cytokines (Reynolds et al., 2006) . The bifurcation analysis of the ODE model demonstrates an array of cytokine dynamical behaviours, from stabilizing at a low concentration, bistable states, oscillations around a high cytokine levels, relapse-remission flare-ups, and convergence to a high cytokine level.

In addition, we extend the ODE model to new stochastic differential equation (SDE) models. The first SDE accounts for local variability in cytokine levels at the initiation of the infection process (cell level), referred to as the demographic variability; whereas the second SDE model is at a larger spatial scale (tissue level), referred to as the environmental variability.

In Section 2, we first summarize the ODE Baker model by interpreting the parameters in terms of cytokine activation and four cytokine regulation mechanisms. We then explain the Hill coefficients as a cooperative binding in cytokine secretion and introduce the two new SDE models. In Section 3, a bifurcation analysis is carried out to identify the bistable region, in which relapse-remission or flare-up patterns are generated. The dynamics of the models are illustrated through both symbolic computation and numerical simulations. In Section 4, the analytical and numerical results are interpreted in terms of the immune response typical of natural reservoirs or of human hosts.

A mechanistic approach to modelling gene expression and/or protein dynamics often begins with an interaction graph to capture pairwise interactions between components and then develops them into a Boolean network model. These models can then be converted to a system of ODEs to help reveal dynamical behaviour of each component (Samaga & Klamt, 2013; Shao et al., 2015) . Hill functions are a common way to transform Boolean functions to ODEs (Hill, 1910; Wittmann et al., 2009; Samaga & Klamt, 2013) .

Here, we apply the Baker model which applies Hill functions to capture the rate of the interactions occurring at the molecular level. Hill functions are often used to model cytokine production (Chaudhry et al., 2004) , which have the general functional form

These functions model the rate of binding of a macromolecule to a ligand, expressed in terms of ligand concentration x. The Hill coefficient m represents the degree of cooperative binding, whereas the parameter k denotes ligand concentraton where half of the binding sites are occupied. In the Baker model, described below, the Hill functions are given in (2.2).

In addition to the Baker model, several other models that have been applied to the study of cytokine dynamics include Seymour & Henderson (2001 ), Yiu et al. (2012 , Jarrett et al. (2014) , Waito et al. (2016) and Waters et al. (2018) . Seymour & Henderson (2001) modelled cytokine concentrations and regulation mediated by cell surface receptors for two specific pro-inflammatory and anti-inflammatory cytokines, IL-1 and IL-10. Yiu et al. (2012) fit a linear ODE model to patient data on serum levels of several pro-inflammatory cytokines, representative of a cytokine storm. Jarrett et al. (2014) developed a 4D model for a bacterial-induced pro-and anti-inflammatory response. The model captured different experimental outcomes. In another ODE model, the unbalanced interaction between several specific proand anti-inflammatory cytokines that induce a cytokine storm is demonstrated by Waito et al. (2016) . Waters et al. (2018) focused their ODE model on IL-2 and its receptor IL-2R in the immune response.

According to their functional roles, we categorize cytokine compounds as pro-inflammatory or anti-inflammatory, denoted by p and a in (2.1), respectively. The generation of both pro-and antiinflammatory cytokines is stimulated by pro-inflammatory cytokines and shows a positive feedback regulation. This regulation is captured by the Hill functions in (2.2). The positive feedback between these two cytokines is denoted as φ(p) and ψ(p) in (2.1), The inhibition effect of anti-acting on proinflammatory production shows a negative feedback and is modelled as θ(a) in (2.1).

The Baker model provides a simple framework for describing the cytokine dynamics. Their model is described in more detail in Section 2 in Baker et al. (2012) . The differential equations are dp dt

where c i > 0 and m i is a positive integer, i = 1, 2, 3. More specifically, the values of m i in (2.1) are determined by cooperative binding to receptors of secreted cytokines. Parameter c 3 stands for the efficiency of this inhibition. Parameters c 1 and c 5 denote the maximum production rates for pro-and antiinflammatory cytokines due to pro-inflammatory binding, while c 0 is the background pro-inflammatory production rate. Parameters c 2 , c 4 and c 6 represent the corresponding cytokine concentrations at which the Hill functions φ( p), θ(a) and ψ( p) reach half of their maximum value, respectively (Santillán, 2008; Stefan & Le Novère, 2013) . Non-dimensionalizing the Baker model (2.1) gives the dimensionless variables p * , a * and t * as p = p * c 2 , a = a * c 4 and t = t * /d a . For notational simplicity, the asterisks are omitted, and the system (2.1) reduces to a model with only five parameters, γ and α i , i = 1, 2, 3, 4. dp dt

The five dimensionless parameters, expressed in terms of the original parameters, are

(2.4)

The new dimensionless parameters α 1 , α 2 , α 4 and γ are scaled by the clearance rate of anti-inflammatory cytokines. Parameter γ is the ratio of clearance rates for pro-inflammatory to anti-inflammatory. In general, increases in α 1 and α 2 upregulate the pro-inflammatory response, whereas an increase in α 4 upregulates the anti-inflammatory response.

Cytokine levels are governed by several mechanisms. Four mechanisms are summarized by Fernandez-Botran (1991) as regulation of (i) cytokine secretion, (ii) cytokine receptor expression, (iii) modulation of one cytokine by other cytokines and (iv) soluble cytokine-binding factors and/or inhibitors. Mechanism (i) affects the pro-and anti-inflammatory cytokines secretion and alters the parameters c 1 and c 5 or α 2 and α 4 in (2.4). The receptor expression in mechanism (ii) is modelled as Hill coefficients in (2.1) as m i , i = 1, 2, 3. The regulation strength of anti-inflammatory cytokines acting on pro-inflammatory cytokines in mechanism (iii) is c 3 or α 2 in (2.4). Binding of soluble cytokine receptors to membrane receptors may downregulate pro-inflammatory cytokines and increase d p or γ in (2.4). Soluble IL-6 (sIL-6) binds with IL-6 and promotes binding to the receptor. This results in an upregulation and increase in the value of parameter c 0 or α 1 in (2.4).

The four dimensionless parameters α 1 , α 2 , α 4 and γ in (2.4) are the focus of the bifurcation analysis. Parameter α 3 is the assumed constant, α 3 = 0.5. As shown in the analysis of the Baker model, variations in the parameter α 3 exhibit similar bifurcation dynamics as γ and α 1 . Variation of the remaining parameters is interpreted in terms of the four mechanisms. Since elevated pro-inflammatory cytokine concentrations are driven by the parameter α 2 , this parameter is chosen as the bifurcation parameter. The responses elicited by changes in the dimensionless parameters are summarized below: α 1 : pro-inflammatory upregulation by soluble cytokine receptors and other background factors, its increase promotes the inflammatory response;

α 2 : pro-inflammatory self-upregulation, its increase promotes the inflammatory response;

α 4 : anti-inflammatory upregulation, its increase inhibits the inflammatory response;

γ : ratio of the pro-inflammatory to anti-inflammatory clearance rate, extracellular and soluble cytokine receptors may increase the ratio, its increase in turn down regulates the inflammatory response.

When cytokine proteins bind to cell surface receptors, this upregulates intracellular signalling and alters the cell function. Both autocrine and paracrine signalling may lead to upregulation or downregulation of pro-and anti-inflammatory cytokines as modelled through the functions φ(p) and ψ(p) in (2.1). Likewise, the upregulation of pro-inflammatory cytokines by anti-inflammatory cytokines is modelled as θ(a) in (2.1). The values of the Hill coefficients m i , i = 1, 2, 3 in (2.1) model the amount of cooperative binding. Instead of representing the number of binding sites, Santillán (2008) states 'that the Hill coefficient is more appropriately described as an interaction coefficient, reflecting cooperativity'. If m i > 1, there is a positive cooperativity, while if m i < 1, there is a negative cooperativity (Stefan & Le Novère, 2013) . The binding of pro-inflammatory cytokines to cell surface receptors results in an increase of the production of both types of cytokines and may also increase the cytokine receptors on the surface of the binding cell. Therefore, we choose m i = 2, i = 1, 2, 3 in (2.1) to model the cooperativity and the sigmoid shape of the regulatory interactions of the Hill functions.

The ODE model in (2.1) serves as a framework for formulation of stochastic models that account for variability in local cytokine activity at the cellular level and the activity at a large scale such as the tissue level. The first SDE model represents the early stage of infection. At this stage, the infection is at the cellular level which is modelled assuming variability in production from individual cells, a demographic variability. With the progression of the infection, the second SDE model represents the infection at the tissue level, a larger spatial scale, assuming environmental variability. The assumption is that pro-inflammatory and anti-inflammatory cytokine production rates, α 2 and α 4 , are affected by factors external to the local environment as the result of recruitment of other cells. Both SDE models apply to either the natural reservoir or the human host. It is the specific parameter region that determines whether the model applies to the natural reservoir or the human host. 

cytokines, we derive a system of SDEs from first principles. System (2.1) is used as a framework for formulation of a birth (upregulation) and death (downregulation) process. The units of p and a are on the order of picomoles per milliliter (Yiu et al., 2012) . Let V denote a fixed volume in which these reactions take place. For example, in an in vitro study, the volume may be on the order of milliliters. We also need to convert picomoles to molecules using Avogadro's constant 6.022 × 10 23 , which is the number of molecules per mole. Therefore, the number of molecules per picomole is ζ = 6.022 × 10 11 . To put the units of p and a in terms of molecules, multiply p and a by ζ V in model (2.1). Let X = (p, a) T be a vector of continuous random variables and Δt > 0 be sufficiently small so that at most one birth or one death occurs during the time interval Δt. Let Δ X = X(t + Δt) − X(t). The possible changes in Δ X are summarized in Table 1 . Note that Δ X(t) depends on the number of molecules in time Δt.

To order Δt, the expectation equals

and the covariance matrix is

Following a similar argument as in Allen (2007) and Allen et al. (2008) , we have a system of Itô SDEs:

where W 1 (t) and W 2 (t) are two independent Wiener processes. To express the SDE model in the same dimensionless variables as in (2.1), first, we apply the same scaling as for the ODE model, p = p * c 2 , a = a * c 4 and t = t * /d a . Note also that dW(t) = √ dt η, where η ∼ N(0, 1). Therefore,

With the preceding substitutions, the units of p * and a * agree with the ODE model (2.1). Second, to remove the additional scaling ζ V introduced in the formulation of the SDE model, we divide the equations by ζ V and drop the * notation for simplicity. The new SDE model has the following form:

with Hill coefficients m i = 2 for i = 1, 2, 3. The additional factors of 1/ √ ζ V that appear in the SDE model (2.5) account for the volume in which interactions take place. If the value of ζ V is sufficiently large, the dynamics are nonlocal and the SDE model (2.5) reduces to the ODE model (2.1). The spatial range of cytokine communication may be on the order of micrometers. The range depends on several factors, including the cell density, the receptor number, the rate of cytokine production and on whether the signalling is autocrine or paracrine (Thurley et al., 2015) . As model (2.5) is a diffusion approximation of a Markov chain model (Table 1) , solutions of (2.5) may become negative. However, numerical solutions are restricted to p 0 and a 0, by setting the numerical values equal to zero for those that reach a negative value.

In the second SDE model, we assume that the production rates of the cytokines vary in time. The variability may be due to heterogeneity in individual cell responses, the local environmental responses formed through the stochastic contributions of each cell (infected, standby and infiltrating immune cells). The parameters α 2 ≡ α 2 (t) and α 4 ≡ α 4 (t). These two parameters promote or inhibit the pro-inflammatory response, respectively. To model the environmental variability, we use SDEs which follow a mean reverting process (Allen, 2016; Iacus, 2009) ,

(2.6)

The constant r i > 0 is the rate of return to the mean rateᾱ i through the body's natural regulation or through therapy (i = 2, 4). During the return process, deviations from the mean are modelled through a term proportional to α i , i.e. σ i α i (t), where σ i > 0. Solutions to (2.6) are non-negative (Iacus, 2009 ). There are other forms for a mean-reverting process that could be considered that also have the property of being non-negative, such as the Cox-Ingersoll-Ross model, where the term with σ i is replaced by

The second SDE model consists of the two SDEs in (2.6) coupled with the differential equations for p and a in model (2.3). The differential equations in (2.3) and in (2.6) form a system of SDEs. The first and second moments of α i can be easily computed by applying Itô's formula (Allen, 2010) . Let

The asymptotic behaviour of the mean and variance of α i in the case of rapid return rate r i is

As the rate r i → ∞ in (2.7), the return to a constant level,ᾱ i , through individual response or via therapy is extremely rapid. In the limit, the constant cytokine production rateᾱ i is maintained. The asymptotic behaviour as t → ∞ of the mean and variance of α i is

(2.8)

If the return rate satisfies r i > σ 2 i /2, the mean production rate approaches a steady-state level with a constant variance as shown in (2.8). (An example of the limiting density for α 2 (t) is graphed in the appendix for an example satisfying σ 2 = √ r 2 .) If the return rate is too slow, r i σ 2 i /2, then the cytokine production rate α i will have large variability. In this latter case, the production rate cannot be maintained at a steady-state levelᾱ i .

Cytokine dynamics are controlled by many regulatory mechanisms. We consider the regulation from anti-inflammatory cytokines and soluble cytokine receptors. Increase in the values of α 4 and γ has an inhibitory effect on pro-inflammatory cytokines. In contrast, the pro-inflammatory response is upregulated with an increase of the rates α 2 and α 1 . This upregulation can occur via soluble cytokine receptors in the small neighbourhood surrounding the cell that facilitate cytokine binding. Since elevation of pro-inflammatory cytokines and dysregulation of the elevated cytokines returning to homostasis are primary causes of a cytokine storm, the parameter α 2 is considered as the bifurcation parameter for the remainder of the analysis. In addition, we consider γ , α 4 and α 1 as control parameters in the following subsections.

Setting γ as the control parameter in Section 3.1, we obtain four regions in γ -space that illustrate the effect of inhibition from γ . These four γ regions are further subdivided into regions where either oscillations are present or absent. Then, the γ value is fixed in each of the two oscillation regions and α 4 is designated as a control parameter. In these two cases, we investigate the effect of the inhibition from anti-inflammatory cytokines. The results are described in Sections 3.2 and 3.3. Finally, in Section 3.4, under intense inhibition from both γ and α 4 , the parameter α 1 is considered as a control parameter to reveal relapse-remission cytokine dynamical behaviour.

As in the Baker model, we assume that there is cooperative binding with m i = 2, i = 1, 2, 3. We verify the fact that if m 1 > 1, model (2.1) can undergo a Hopf bifurcation, which potentially induces oscillations and even relapse-remission behaviour with high level spikes in pro-inflammatory concentration. Such type of behaviour occurs in recurrent infections (Zhang et al., 2014b) and relapseremission autoimmune diseases (Zhang et al., 2014a) .

The bifurcation analysis delimits parameter regions as boundaries. Close to these boundaries, small changes in parameter input can cause a dramatic effect in the model output. Moreover, this bifurcation We illustrate the dynamics of the SDE models when variability is included at the cell or tissue level. To illustrate the local variability at the cell level, in the SDE model (2.5) up to time t = 50 and we let ζ V = 1. For the widespread variability at the tissue level, in SDE (2.6), we assume that the mean of the activation ratesᾱ 2 andᾱ 4 is the same as in the ODE model (2.3), α 2 and α 4 . In addition, we assume that σ i = √ r i , i = 2, 4. Thus, it follows from (2.8) that the mean and asymptotic variances satisfy

The preceding assumption implies that the limiting mean and variance do not depend on r i . The rate of return is assumed to be r i = 0.5, i = 2, 4, for the simulations. In this subsection, we fix parameter α 4 = 3.5 at a low value, while γ and α 2 are free. The 2D bifurcation dynamics in α 2 -γ space are summarized in Fig. 1(a) . The cytokine dynamical behaviour is demonstrated for three γ values from non-Hopf regions (i) and (ii), two-Hopf region (iii) and one-Hopf region (iv).

In the non-Hopf region (i) and (ii), fixing γ = 1.25, Fig. 1(b) shows that model (2.3) experiences two SN bifurcations at SN 1 and SN 2 . Fixing γ = 1.28 in the two-Hopf region (iii), Fig. 1(c) shows that model (2.3) exhibits oscillation when α 2 ∈ (Hopf 1 , Hopf 2 ). In addition, model (2.3) exhibits bistability between two SN bifurcations, SN 1 and SN 2 . In the bistable interval, there exists a pro-inflammatory oscillation. Further increases in γ , such as γ = 2 as shown in Fig. 1(d) , increase the bistable region between the two SN bifurcations. As the value of α 2 increases, the pro-inflammatory cytokines first stabilize at a healthy low level on the left of SN 2 , then to a bistable state with α 2 ∈ (SN 2 ,SN 1 ). In the bistable region, under outside stimulation the pro-inflammatory cytokines have the possibility to return to a healthy low level. However, the possibility of returning to the healthy state disappears when α 2 increases and passes SN 1 .

Taking three values of γ = 1.25, 1.28, 2, the corresponding 1D bifurcation diagrams are shown in Fig. 1(b-d) . Then we fix α 2 = 20, e.g. the value of the equilibriump is downregulated slightly with increasing of the γ value. The strong pro-inflammatory response can be reached with a large value of α 2 in all three cases. However, the pro-inflammatory response may vary periodically. This phenomenon is shown as bistability with oscillation from a Hopf bifurcation, in Fig. 1(d) , where α 2 ∈ (SN 2 , Hopf). For the SDE models, the bistable region may be a representative of the inflammatory response in a natural reservoir, whereas the region of high pro-inflammatory response may represent the inflammatory response in a human infection.

For the bistable region, consider parameter values corresponding to Fig. 1(b) . For α 2 = 10, α 4 = 3.5 and γ = 1.25, the ODE exhibits bistability. However, in Fig. 2 , the variability in the two SDE models shows that the sample paths fluctuate between the two stable states. For the SDE model (2.5), there are large fluctuations in both pro-and anti-inflammatory cytokines. However, for SDE (2.6) the response is more controlled; the cytokine levels are low but persistent with occasional flare-ups, such as in the infection of a natural reservoir. The pro-inflammatory cytokines are slightly upregulated through variations in the parameter α 2 but the regulatory response is weak due to low values of α 4 and γ .

A strong pro-inflammatory upregulation with a large α 2 = 90 and weak regulation through low α 4 = 3.5 are seen Fig. 3 . The ODE model is in the region of a high stable pro-inflammatory response ( Fig. 1(d) ). At very low cytokine levels, the ODE model overshoots the equilibrium value before returning to the stable equilibrium, Fig. 3(a) . Thus, in the two SDE models, when variability reduces the dynamics to low levels, a return to equilibrium may also result in an overshoot phenomenon. In the two SDE models, there is a variability at both the cell and tissue levels, driven by the large proinflammatory response through α 2 but a weak anti-inflammatory response through α 4 . Such behaviour may be indicative of the cytokine storm seen in human infection. In the appendix, the dynamics of SDE model (2.6) are explored for different values of the rate of return r = r i , i = 2, 4 to the mean value of α i , i = 2, 4. For short time spans and small values of r, the α i values change slowly; but for large values, the change is more rapid.

3.2 Parameter α 4 vs. α 2 for small γ In this subsection, inhibition from anti-inflammatory cytokines is investigated under a small value of γ . The value of γ = 1.25 is taken in region (ii) in Fig. 1(a) . The generation rate of anti-inflammatory Fig. 3 . Dynamics of the ODE and SDE models with ODE parameter values equal to α 2 = 90, α 4 = 3.5 and γ = 2. The ODE model is in the region of a high stable pro-inflammatory level, α 2 > Hopf, α 2 ( Fig. 1(d) ). The ODE model solution overshoots the equilibrium when initial conditions are small, p(0) = 0.1, a(0) = 0.1. The two SDE models exhibit large variability in cytokine levels about the equilibrium, indicative of a cytokine storm. cytokine, α 4 , is chosen as the control parameter with bifurcation parameter α 2 . Other parameter values are fixed as in Table 2 . The bifurcation analysis is illustrated in α 2 -α 4 parameter space in Fig. 4 .

The bistable region, shown in the dotted region of Fig. 4(a) , narrows down and diminishes with the increase of the α 4 value. Nevertheless, taking α 2 = 20, e.g. Fig. 4(b-d) shows that the large value of α 4 decreases the value of the equilibriump dramatically. Therefore, with a weak downregulation such as from soluble cytokine receptors, anti-inflammatory regulation reduces the pro-inflammatory response to a low healthy state.

Similar behaviour for the ODE and SDE models, as exhibited in Fig. 2 , holds for parameter values chosen from the bistable region, α 2 ∈ (SN 2 , SN 1 ) in Fig. 4(b) . The regions where there is a greater anti-inflammatory response, large α 4 exhibit periodicity and are explored further in the next section.

Regulation from anti-inflammatory cytokines is studied with a large value of γ , i.e. in region (iv) of Fig. 1(a) . Specifically, we take γ = 1.5 and α 4 is chosen as the control parameter with α 2 the bifurcation parameter. See Fig. 5 .

Bifurcation curves in Fig. 5 (a) delimit the α 2 -α 4 parameter space into three parts: (i), (ii) and (iii) in Fig. 5(a) . In region (i), the value of α 2 is small, pro-inflammatory response converges to a low-level healthy state with a moderate level of α 2 and intensive inhibition regulations from both antiinflammatory cytokines and soluble cytokine receptors. In contrast, region (iii) implies a cytokine storms with the amplifying α 2 values. The region (ii) in Fig. 5(a) is an intermediate stage between (i) and (iii), where the pro-inflammatory response oscillates between low and high levels. The oscillations spend more time at the low pro-inflammatory level with larger values of anti-inflammatory response α 4 . Under a moderate pro-inflammatory generation rate α 2 , which locates close to and on the right side of SN 1 in Fig. 5(c) , the oscillating behaviour shows as a relapse-remission pattern of pro-inflammatory response. In the α 2 -α 4 parameter plane, model (2.3) shows bistability and oscillation in the dotted and grey regions, respectively. Furthermore, taking α 4 = 3, 8, 28 individually, F 2 (p, α 2 ; α 4 ) = 0 is plotted as 1D bifurcation diagram in (b), (c) and (d). Solid and dashed curves denote stable and unstable equilibria. The unstable equilibrium branch between SN 1 and SN 2 separates the initial condition space for bistability, while the unstable branch between Hopf 1 and Hopf 2 is enclosed by stable periodic solutions.

Further increases in the anti-inflammatory response α 4 and in the clearance rate of pro-inflammatory cytokines through γ , the ODE model exhibits oscillatory dynamics (Fig. 6 with α 2 = 20) . The oscillations are not evident in the two SDE models. The existence of the stable equilibrium at a low pro-inflammatory response allows the cytokine dynamics in the SDE model (2.6) to return to the low state for an extended period of time. The dynamics are similar to those in Fig. 2 (infection in a natural  reservoir) .

For a larger pro-inflammatory response α 2 = 50, the anti-inflammatory cytokines cannot control the pro-inflammatory response, resulting in a cytokine storm. The dynamics in Fig. 7 are similar to those in Fig. 3 . The parameter α 1 represents background production of pro-inflammatory cytokines. Strong cytokine responses may influence this background pro-inflammatory cytokine production. For example, soluble cytokine receptors facilitate cytokine signalling, such as α2-macroglobulin (James, 1990 ) carry IL-1 and IL-6 in circulation (Fernandez-Botran, 1991) and increase the background pro-inflammatory cytokine concentration. We consider parameter α 1 as a control parameter, with a bifurcation analysis in the α 1 -α 2 parameter space, a parameter region not fully explored in the Baker paper. See Fig. 8 .

Numerical simulations show that relapse-remission behaviour occurs when the parameter pair (α 1 , α 2 ) locates close to the SN curve on its right side with a moderate value of α 2 ∈ (20, 40). Thus, relapse-remission pattern in pro-inflammatory response is more likely to occur when the background pro-inflammatory production rate, α 1 , increases under a moderate pro-inflammatory production rate, α 2 . There is an even more vigorous anti-inflammatory response with high values of γ and α 4 . The dynamics of the ODE model are illustrated for two different sets of parameter values in Fig. 8(b and c) . The corresponding stochastic dynamics are illustrated for one set of parameter values in Fig. 9 .

The pro-inflammatory response is not completely controlled due to this background proinflammatory production. Therefore, moderate pro-inflammatory production rate α 2 can drag the pro-inflammatory concentration to high levels periodically. This interaction between promoting and inhibiting regulations forms spikes in pro-inflammatory concentrations, which lead to high viral replication and release illustrates a persistent infection with flare-ups seen in a natural reservoir. Fig. 8(b) . The high antiinflammatory response is able to control the pro-inflammatory response for short periods of times, but the underlying oscillations show occasional flare-ups, indicative of the natural reservoir.

The goal in this investigation is to apply simple mathematical models for pro-and anti-inflammatory cytokines to describe the dynamics of viral zoonotic infectious diseases in different hosts. Through bifurcation analysis and numerical simulations, the models' outcomes are related to the responses seen in natural reservoirs and in human hosts. In natural reservoirs, the infection persists and the inflammatory response is controlled by immune homeostasis through largely undefined mechanisms. In human hosts, the infection can result in a stronger dysregulated pro-inflammatory response which eventually leads to tissue damage.

The regulatory interaction between pro-and anti-inflammatory cytokines is captured by a 2D ODE system. The positive feedback in cytokine generation and the negative feedback inhibition where antiacts on pro-inflammatory cytokines are modelled through Hill functions. The positive cooperativity in cytokine bindings and the non-linear sigmoid shape of the regulatory interaction are modelled by the Hill coefficient, whose magnitude is greater than unity. The underlying ODE model exhibits a broad range of behaviours from low to high-stable cytokine levels, bistable states, oscillations and relapse-remission flare-ups. Regions in parameter space are associated with upregulation of either pro-and/or anti-inflammatory cytokines. The bifurcation analysis is carried out through symbolic computation. This parameter-free approach yields algebraic formulas for use in qualitative analysis. The bifurcation diagrams in Figs 1, 4, 5 and 8 reveal parameter regions where there exists a single stable equilibrium or bistability or periodic solutions in γ − α 2 , or α 2 − α 4 , and α 1 − α 2 space. For example, the low pro-inflammatory generation rate, α 2 , may keep the pro-inflammatory cytokines level, p, at a low stable steady state. This is illustrated by the lower branch of p on the left of SN 2 in Fig. 1(b-d) . While the strong pro-inflammatory cytokine responses, a high level of p occurs on the right of SN 1 in Fig. 1(b-d) . The bistable region, between SN 1 and SN 2 , increases with the growth of the decay ratio γ between pro-and anti-inflammatory cytokines.

As the ODE model does not capture the variability at the cell and tissue levels, the ODE model is extended to two SDE models. The numerical examples illustrate that some of the distinct outcomes associated with the ODE model cannot be distinguished in the SDE models. For example, bistability or oscillatory behaviour in the ODE is seen in the SDE models as a chronic infection with occasional flareups, typical of infection in the natural reservoir (Figs 2, 6 and 9) . Also, the high-stable cytokine levels may be viewed as a cytokine storm in the SDE models (Figs 3 and 7) . The flare-ups illustrated in both the ODE and SDE models are also characteristic of relapse-remission immune responses seen in chronic infection. This type of behaviour has been demonstrated in other types of ODE models (Zhang et al., 2014a (Zhang et al., ,b, 2016 . The examples in Sections 3.1-3.4 show the control of the pro-inflammatory response through large γ or large α 4 , resulting in low levels of p and a but occasional flare-ups (bistable or oscillatory regions), indicative of the natural reservoir. The examples also illustrate cases, where the proinflammatory response α 2 cannot be reduced to low levels, through γ or α 4 (high endemic equilibrium for p), indicative of a cytokine storm in human infection.

In summary, we proposed a simple two-component model for a viral infection of cells or organs (e.g. lung, liver) to gain insight into the effects of two opposing inflammatory responses that give rise to a cytokine storm. The theoretical mathematical models are preliminary and non-mechanistic. They illustrate potential outcomes of pro-versus anti-inflammatory cytokines through four model parameters with broadly defined functional roles (described in Section 2.2, (i)-(iv)).

However, the mechanisms responsible for the cytokine storm are much more complex. To develop a mechanistic model requires time series data from carefully designed laboratory experiments that identify when specific pathways and genes/proteins are induced or suppressed. Hence, these simple models provide a framework for such experimental data. In other words, these complex data can be used to test some outcomes of these models. The values for the parameters andp should satisfy a necessary condition as follows:

F h = − α 2 4 + 1 (γ + 1)p 9 − 2 α 2 3 (γ + 1) + α 2 4 + 1 p 7 − 2 α 1p 6 + −α 4 3 + α 2 4 + 1 γ − α 2 3 α 2 3 + 4 − α 2 4 − 1 p 5 − 2 2 α 2 3 + 1 α 1p 4 − 2 α 2 3 α 2 3 − γ + 1 p 3 − 2 α 2 3 + 2 α 1 α 2 Fig. A1 . The graph of limiting pdf of the random variable α 2 (t) in the SDE model (2.6) with mean equal to 90. Ultimately, the variable α 2 (t) takes on the full range of values associated with the pdf. sample paths for p(t) and a(t) as illustrated in Fig. A2 . The rate of change for case (c) r = 5 is much faster than for case (a) r = 0.05. As a final note, the numerical simulations are for the dimensionless ODE and SDE models. Data on the actual response time and magnitude of the immune response for individual hosts will help in parameterization of models to differentiate natural reservoirs versus human hosts.

Modeling with Itô Stochastic Differential Equations

Environmental variability and mean-reverting processes

An Introduction to Stochastic Processes with Applications to Biology

Construction of equivalent stochastic differential equation models

Mathematical modelling of cytokine-mediated inflammation in rheumatoid arthritis

Proinflammatory cytokine responses induced by influenza a (h5n1) viruses in primary human alveolar and bronchial epithelial cells

Empirical models of the proliferative response of cytokine-dependent hematopoietic cell lines

Seoul virus enhances regulatory and reduces proinflammatory responses in male Norway rats

Regulatory T cells enhance persistence of the zoonotic pathogen Seoul virus in its reservoir host

Soluble cytokine receptors: their role in immunoregulation

SARS coronavirus and innate immunity

Expression of elevated levels of proinflammatory cytokines in SARS-CoV-infected ACE2+ cells in sars patients: relation to the acute lung injury and pathogenesis of SARS

Soluble cytokine receptors

The cytokine network in infectious diseases

The possible effects of the aggregation of the molecules of haemoglobin on its dissociation curves

Bat origin of human coronaviruses

Interactions between cytokines

Modelling the interaction between the host immune response, bacterial dynamics and inflammatory damage in comparison with immunomodulation and vaccination experiments

Global trends in emerging infectious diseases

A global perspective on hantavirus ecology, epidemiology, and disease

Bats are natural reservoirs of sars-like coronaviruses

The cytokine storm of severe influenza and development of immunomodulatory therapy

Epidemic dynamics at the human-animal interface

Role of vascular and lymphatic endothelial cells in hantavirus pulmonary syndrome suggests targeted therapeutic approaches

Cytokines in the pathogenesis of rheumatoid arthritis

Persistent hantavirus infections: characteristics and mechanisms

Innate immune responses to influenza A H5N1: friend or foe?

Innate immune response of human alveolar type ii cells infected with severe acute respiratory syndrome-coronavirus

A reduced mathematical model of the acute inflammatory response: I. derivation of model and analysis of antiinflammation

Pathogenesis and host response in Syrian hamsters following intranasal infection with Andes virus

Modeling approaches for qualitative and semi-quantitative analysis of cellular signaling networks

On the use of the hill functions in mathematical models of gene regulatory networks

Kinetics of immune responses in deer mice experimentally infected with Sin Nombre virus

Pro-inflammatory-anti-inflammatory cytokine dynamics mediated by cytokine-receptor dynamics in monocytes

From boolean network model to continuous model helps in design of functional circuits

A comparative analysis of host responses to avian influenza infection in ducks and chickens highlights a role for the interferon-induced transmembrane proteins in viral resistance

Cooperative binding

Role of host cytokine responses in the pathogenesis of avian H5N1 influenza viruses in mice

Three-dimensional gradients of cytokine signaling between t cells

Into the eye of the cytokine storm. Microbiol

Zika virus, vectors, reservoirs, amplifying hosts, and their potential to spread worldwide: what we know and what we should investigate urgently

A Mathematical Model of Cytokine Dynamics During a Cytokine Storm

Emerging zoonotic viral diseases

Avian influenza viruses, inflammation, and CD8+ T cell immunity

The effects of interleukin-2 on immune response regulation

Influenza A viruses target type II pneumocytes in the human lung

Transforming Boolean models to continuous models: methodology and application to T-cell receptor signaling

Bats and viruses: friend or foe?

Modeling and analysis of recurrent autoimmune disease

Viral blips may not need a trigger: how transient viremia can arise in deterministic in-host models

Backward bifurcations, turning points and rich dynamics in simple disease models

Herep satisfies F 4 p; α 2 , α 1 = 0 in (A.14). A BT bifurcation occurs at α 1 , α 2 = (0.0286, 15.4864). In the 2D bifurcation diagram in Fig

A.5 Probability densities for parameters in the SDE model

In this case, the shape of the limiting pdfs only depends on the mean value. See Fig. A1. Sample paths of the random variables p(t) and a(t) for the SDE model (2.6) are graphed in Fig. A2 for different values of r = r 2 = r 4 (and other parameter values as in Fig. 3). The return rate r differs in the three graphs below, r = 0.05, 0.5 and 5. The limiting mean value for the parameter α 2 (t) is equal to 90 and for α 4 (t) it is 3.5. The rate of return impacts the rate of change of the parameters α i (t), i = 2, 4

We thank the referees for their helpful comments and suggestions on the original manuscript. 

We expand on the analysis presented for the Baker model. Setting the left-hand side of (2.3) to 0, we obtain equilibrium solutions (p,ā) of model (2.3), wherēandp satisfies the following equation:The bifurcation dynamics depend on the local stability of the equilibria, (p,ā), since stable (unstable) equilibria attract (repel) nearby solutions. The corresponding eigenvalues of the Jacobian matrix of system (2.3) determine the stability of equilibria and are obtained from the roots of the characteristic polynomial λ 2 − Trace (J) λ + Det (J) = 0, whereThe stability of an equilibrium changes when the system undergoes a bifurcation. An SN bifurcation occurs when Det . For a dynamical system undergoing bifurcation, a small perturbation or uncertainty in parameter values can cause a dramatic change in model behaviour. Therefore, bifurcation curves in terms of parameter values form the boundaries for different model behaviours.The value of m i , i = 1, 2, 3 is related to the number of binding sites on the surface of cytokinesecreting cells. Generally, larger values of m i give rise to richer dynamical behaviour of models (2.1, 2.3). Two static bifurcations usually give rise to bistability. Hopf bifurcation plays an important role as a source for oscillations. Therefore, we investigate the necessary condition for a Hopf bifurcation. Considering F (p) = 0 of (A.1) into Trace (J) = 0 in (A.2), we have. Therefore, Hill coefficients are set as the smallest integers which give rise to rich dynamical behaviours that is m i = 2, i = 1, 2, 3. Moreover, since Hopf bifurcation is a codimension-one bifurcation, the dynamical behaviour of the model can fully unfold with the variance of one-parameter value. In this study, we concentrate on investigating how the extra pro-inflammatory cytokinesecreting rate, α 2 , causes cytokine storm in virus infections. Therefore, α 2 is chosen as a bifurcation parameter. The bifurcation curves (the corresponding equilibrium undergoes one specific bifurcation when parameter values locate on one of these curves) are determined by The bifurcation analysis shows that an SN bifurcation occurs if Det 1 p; γ , α 2 = 0 in (A.7); Hopf bifurcation occurs if Trace 1 p; γ , α 2 = 0 in (A.6) and Det 1 p; γ , α 2 > 0. In (A.7) and (A.6),p satisfies equation F 1 p; γ , α 2 = 0 in (A.5). (A.7) and (A.6) give the SN and Hopf bifurcation curves, denoted as SN 1 , SN 2 and Hopf in Fig. 1 . These bifurcation curves divide the α 2 -γ parameter space in Fig. 1 Here, γ hm = 1.2765 denotes the minimum γ value for the occurrence of Hopf bifurcation, and γ BT = 1.2845 denotes BT bifurcation critical point. Moreover, according to (A.4), the necessary condition of a Hopf bifurcation for m 1 = 2 is γ > 1. Therefore, in region (i), a Hopf bifurcation does not occur in the entire parameter space; while in region (ii), Hopf bifurcation does not occur in the current parameter set, but there is the possibility in other parameter sets, investigated in the next subsection. Region (iii) is a transition case. Model (2.3) exhibits oscillations in a slim window of α 2 -values, sandwiched between two Hopf bifurcations. The eigenvalues of one of the two Hopf bifurcations become a pair of real numbers with opposite signs in region (iv). There is a large interval for α 2 values where model (2.3) exhibits oscillations. In Fig. 4(a) , an SN bifurcation occurs if Det 2 p; α 2 , α 4 = 0 in (A.10); a Hopf bifurcation occurs if Trace 2 p; α 2 , α 4 = 0 in (A.9) with Det 2 p; α 2 , α 4 > 0 in (A.10). The SN and Hopf bifurcation conditions give the corresponding curves as SN 1 , SN 2 , Hopf 1 and Hopf 2 in Fig. 4. Model (2.3) shows bistability in the dotted region delimited by two bifurcation curves, SN 1 and SN 2 and oscillation in the grey region delimited by two Hopf bifurcation curves, Hopf 1 and Hopf 2 . Further, 1D bifurcation analysis is carried out for three α 4 values taken at 3, 8 and 28. The corresponding bifurcation diagrams are plotted in Fig. 4(b-d) , respectively. The 1D bifurcation analysis shows that with the increase of the value of α 4 , the bistable α 2 interval between SN 1 and SN 2 shrinks down until it disappears in the case of α 2 = 28 in Fig. 4(d) . The middle branch unstable equilibrium is a saddle and serves as a boundary for bistability.In the top branch equilibrium, with an increase in α 4 , a Hopf bifurcation appears; the distance between the two Hopf bifurcations enlarges and the unstable equilibrium is enclosed by a stable periodic solution. If the lower branch equilibrium exists, it is a stable node. 292 W. ZHANG ET AL.A.3 Parameter α 4 vs. α 2 for large γ Setting α 1 = 0.025 and γ = 1.5 yieldsThe Trace (J) = 0 and Det (J) = 0 in (A.2) considering the equilibrium condition a =ā (p) give Bifurcation analysis gives a Hopf bifurcation curve and SN bifurcation curves, denoted as Hopf, SN 1 and SN 2 in Fig. 5(a) , according to Det 3 p; α 2 , α 4 = 0 in (A.13) and Trace 3 p; α 2 , α 4 = 0 in (A.12) with Det 3 p; α 2 , α 4 > 0, respectively. Herep satisfies equation F 3 p; α 2 , α 4 = 0 in (A.11). The points where the Hopf bifurcation curve is tangent to a SN bifurcation curve give rise to two BT bifurcations at α 2 , α 4 = (17.8016, 15.8314) and α 2 , α 4 = (5.7757, 1.5577).A.4 Parameter α 1 vs. α 2 Setting α 4 = 12, and γ = 1.5, we have F 4 p; α 2 , α 1 =p 7 + α 1 + α 2 p 6 − 873 4p The static bifurcation occurs whenHere in (A.17),p and the parameters satisfy the following condition: F s = α 2 4 + 1 2 γp 15 t + 5 α 2 4 + 1 2 α 2 3 γp 13 + 2 α 2 4 + 1 α 1p 12 + 2 2 α 2 3 α 2 4 + 5 α 2 3 + 4 α 2 4 α 2 4 + 1 α 2 3 − α 2 4 + 1 2 γp 11 + 2 3 α 2 3 α 2 4 + 5 α 2 3 + α 2 4 + 1 α 1p 10 + 3 α 2 4 2 α 4 3 + α 2 4 + 10 α 4 3 − 5 + 2 α 2 4 8 α 2 3 − 1 α 2 3 γp 9 + 6 α 2 3 + 1 α 2 4 + 20 α 2 3 + 10 α 1 α 2 3p 8 + 5 α 4 3 − 10 + 2 4 α 2 3 + 1 α 2 4 α 4 3 γp 7 + 2 α 2 3 α 2 4 + 10 α 2 3 + 3 α 2 4 + 10 α 1 α 4 3p 6 + α 4 3 + 2 α 2 4 − 10 α 6 3 γp 5 + 2 5 α 2 3 + α 2 4 + 10 α 1 α 6 3p 4 − 5 α 8 3 γp 3 + 2 α 2 3 + 5 α 1 α 8 3p 2 − α 10 3 γp + 2 α 1 α 10 3 = 0. (A.18)A Hopf bifurcation occurs when α 2h = α 2hā p 2 α 2 3 +p 2 2 , where α 2ha = α 2 4 + 1 γp 7 − α 1p 6 + 2 α 2 3 + α 2 4 + 1 γp 5 − 2 α 2 3 + 1 α 1p 4 + α 2 3 + 2 α 2 3 γp 3 − α 2 3 + 2 α 1 α 2 3p 2 + α 4 3 γp − α 1 α 4 3 . (A.19)