Modelling the linkage between influenza infection and cardiovascular events via thrombosis

There is a heavy burden associated with influenza including all-cause hospitalization as well as severe cardiovascular and cardiorespiratory events. Influenza associated cardiac events have been linked to multiple biological pathways in a human host. To study the contribution of influenza virus infection to cardiovascular thrombotic events, we develop a dynamic model which incorporates some key elements of the host immune response, inflammatory response, and blood coagulation. We formulate these biological systems and integrate them into a cohesive modelling framework to show how blood clotting may be connected to influenza virus infection. With blood clot formation inside an artery resulting from influenza virus infection as the primary outcome of this integrated model, we demonstrate how blood clot severity may depend on circulating prothrombin levels. We also utilize our model to leverage clinical data to inform the threshold level of the inflammatory cytokine TNFα which initiates tissue factor induction and subsequent blood clotting. Our model provides a tool to explore how individual biological components contribute to blood clotting events in the presence of influenza infection, to identify individuals at risk of clotting based on their circulating prothrombin levels, and to guide the development of future vaccines to optimally interact with the immune system.


S1.1 Model development
We model the influenza infection process using a system of ordinary differential equations (ODEs). Model variables and units are given in Figure 1 Legend. Similarly, parameter values and corresponding descriptions are provided in Supplementary Table S3. A schematic of the model structure illustrating the key model components and their relationships are shown in Figure 1A.
Equations (1)-(3) describe the basic infection process consisting of the viral load, uninfected target cells, and infected target cells. Equations (4)-(12) describe the immune and inflammatory response to infection. In particular, Equation (4) represents IFN-I comprising the key portion of the adaptive immune response and Equation (6) representing CD8 + T cells which is a key component of the adaptive immune system. Equations (1) through (12) are adapted from existing models [21], [22].

S1.1.1 Model for influenza infection, immune response, and inflammatory response.
We model the change in virus with Equation (1). The first term ( ) represents the increase of virus from infected cells, ( ) represents the natural death of virions, and ( ) represents depletion of virions due to infecting target cells. Lastly, the term ( ) represents the inactivation of virions due to antibodies.
We describe target epithelial cells with Equation (2). The first term in this equation represents logistic growth with total target cell population + + limited by carrying capacity 0 . The second term ( ′ ) represents infection of target cells due to virus. Resistant target cells lose protection at rate , hence becoming susceptible to infection. Lastly, target cells become resistant in the presence of IFN-I represented by the last term ( ). Equation (3) represents the change in infected target cell population . Virus infects target cells hence the term ( ′ ) reflects the resulting increase in infected cells. Infected cells decay naturally with death rate . The term ( ) represents IFN-activated natural killer (NK) cells neutralizing infected cells [22]. This is important because NK cells limit the early dissemination of the virus by killing infected cells and secreting inflammatory and antiviral cytokines like IFN− to activate macrophages. Lastly, ( ) models the decay of infected cells by CD8 + T cells, which come from the adaptive immune response [22]. CD8 + T cells play a crucial role in anti-viral immunity, but current influenza vaccines do not lead to a robust and long-lasting CD8 + T cell memory response.
Therefore, focusing our model on these cells provides critical insights into how future vaccines should be developed to improve flu-specific T cell responses. Equation (4) models the innate immune response given by IFN-I ( ). Here infected cells promote the innate immune response and IFN-I decays at rate . Equation (5) represents the resistant target cell population and follows directly from Equation (2). CD8 + T cells (Equation 6) are activated by macrophages, whose growth we model with a Hill function (with Hill coefficient ℎ ). CD8 + T cells are depleted by inactivating infected target cells, reflected in the second term ( ), and decay at rate [21]. Equation (7) models the B cell population. B cells exist naturally, which we capture with the recruitment and decay rates. B cell production, reflected in the term ( − ) is also stimulated by macrophages and IL-10 [21]. Equation (8) models the change in IL-10. The first term captures the production of IL-10 from macrophages. IL-10 (L) and Chemokines (K) are produced in a similar way by macrophages. Due to the receptor saturation in macrophages the Michaelis-Menten saturation term is used to model their growths. IL-10 inhibits its growth by itself in a saturation fashion. The second term ( − (1 − ) ) of (8) represents the decay of IL-10 at rate and production from target cells in the lung, which is limited by presence of . Equation (9) models the change in the dead cells, which follows from equation (3). Note that we only include dead and do not consider damaged epithelial cells in this study. Equation (10) models the change in TNF , denoted by , whose production is stimulated by the interaction of macrophages with the presence of viral particles, exogenous TNF , and dead cells. For the details of the TNF production see [21]. Equation (11) models macrophages (M), which are comprised of i) macrophages which persist in epithelial tissue (i.e. tissue resident macrophages) and ii) macrophages that enter the lung from blood via a chemokine gradient (i.e. inflammatory macrophages) [21]. Tissue resident macrophages are essential for initiating the inflammatory response and recruitment of inflammatory cells to the infection site, including NK cells and monocytes, which ultimately differentiate into inflammatory macrophages. Equation (12) models the production of chemokines (K). Macrophages and pro-inflammatory signal 1 stimulate chemokine production. We also consider the effects of IL-10 and its inhibition of chemokine production [21]. IL-10 mediated regulation of the inflammatory process is important for limiting tissue damage associated with inflammation and initiating tissue repair.
Finally, we note explicitly that TNF production is stimulated by the accumulation of proinflammatory signal Σ2, whose strength increases with the accumulation of dead cells, influenza virus, and TNF itself. Similarly, IL-10 and chemokines are stimulated by the pro-inflammatory signal Σ1. The strength of the signal Σ1 increases with the accumulation of dead cells and influenza virus.

S1.1.2 Model for blood coagulation
Linking inflammation to blood coagulation via TF induction is a key feature of this study. In Figure 1A and Figure 1B we depict the key events considered in this study leading to a blood clot. Inflammatory cytokines, such as tumor necrosis factor TNF , lipopolysaccharides (LPS), IL-1, c-reactive protein, and IFN-I induce the production of tissue factor [21]. Thrombin, which is generated during blood coagulation, has a synergistic effect on TF production [10,12, 18, 40, 41, 49 -51]. On the other hand, anti-inflammatory cytokines such as IL-10 regulate the activity of tissue factor [31]. TF will trigger the blood coagulation process and prothrombin will be consumed. Thrombin is produced as a result of the activation and consumption of prothrombin [14,34]. Meanwhile, thrombin could also enhance the activation of thrombin when a sufficient concentration is present [14,34].
We connect the inflammatory response and blood coagulation via TF induction. Here we model TF production by Equation (13). The first term on the right-hand side of Equation (13) accounts for the synergistic effects of thrombin-TNF and induction from IFN-I. The second term ( [ ][ ]) takes into consideration the activation of prothrombin ([II]) and the last term is natural degradation at rate . Here the growth rate is set to 0 if the TNF concentration is below a threshold . When TNF passes this threshold concentration level, an existing atherosclerotic plaque will be disrupted, resulting in the exposure of TF and formation of a blood clot [52]. With TF induced, blood coagulation begins and thrombin generation occurs. Equations (14)-(16) model the prothrombin, thrombin, and antithrombin dynamics, respectively. We adopt a model of the simplified extrinsic pathway describing thrombin generation or blood coagulation cascade [12]. In particular, we adopt a model of the simplified extrinsic pathway due to the size and complexity of the mechanistic physiological models (e.g. [9,11]). The initiation phase consists of consuming existing prothrombin which forms thrombin at rate . In the propagation phase, thrombin generation from prothrombin occurs at rate . Finally, thrombin inhibition is modelled using the rate constant [12]. The parameters , , are switched on and off based on the amount of thrombin present and as a result are piecewise function (see Supplementary Table S1 for the switching conditions) [12]. Finally, we model the size of a blood clot resulting from active thrombin with Equation (17). We assume, in addition to thrombin presence, the growth of a blood clot requires the presence of TF ([III]). Here represents the percentage of an artery that is restricted by a blood clot. To capture qualitative behavior of blood clot formation, we ensure that is bounded above by 1 and also increases in the presence of thrombin and prothrombin [9].
TF induction threshold T III :

S1.2 Parameter estimation and initial conditions.
Here we include the details of the parameterization of model (1) we develop in the main text. We list model parameters informing the mathematical model (1) in Supplementary Table S3 as well as their respective ranges, units, and sources. Equations (1)-(12) are parametrized using estimates from existing influenza infection modeling studies [21,22]. Equations (14)-(16) are parametrized using existing studies modelling blood coagulation [12,23]. We inform the remaining parameters by integrating several additional sources of experimental and clinical data.

Decay of tissue factor
: In an in vitro experiment, human peripheral blood mononuclear cells (PBMC) were isolated and TF was solubilized from these cells [47]. Using the reported TF halflife of 1.3 hours from this in vitro experiment, we estimate the decay of TF to be 12.9 1/day [47].
Induction of TF: We estimate parameters for Equation (13) in model (1). For the generation of tissue factor via TNF and IFN-I we use in vitro experimental data to fit and ℎ [46]. We also inform the parameters 4 and , which relate the synergistic effects of TNF and thrombin on TF induction, using in vitro data [46]. In particular, these parameters were calibrated to reflect a 6.2-fold increase in TF induction due to TNF when thrombin is present [46]. We note that the data collected from [27,47] was based on healthy subjects. We did not explicitly consider the health status of individuals in this study; but, it is an important factor that may be critical to explore as discussed in the main text as in Section 3.1 as a future direction.
Baseline maximal uptake of TF and substrate affinity ℎ : We proceed using in vitro experimental data measuring TF expression on the surface of epithelial cells [27]. Human umbilical vein cells (HUVECs) were prepared and TF expression was assessed by incubating with TNF and vascular endothelial growth factor (VEGF) for 6 hours [27]. To inform ℎ and we use Equation (13) and the following assumptions: i) no TF decay over the course of the experiment and ii) TNF is constant over the course of the experiment. Also, IFN-is not present in this experiment, so = 0. Hence, we have the following model equation Separating variables results in the following linear equation with variables and ℎ . Now, we use the following data points from the outlined experiment data for TF induced by TNF . We denote the pairs in terms of (TNF , TF) for two experiments: Experiment 1 (0.015 nM, 700 nM) and experiment 2 (0.2 nM, 1450 nM).

Maximal uptake of TF,
: The maximal uptake of TF is enhanced by the presence of thrombin [42]. Here we account for the synergistic effects of thrombin and TNF on TF induction. We use in vitro experimental data which assesses TF expression in HUVECs after incubation with 1) TNF or 2) both TNF and thrombin. Hence, we use the baseline value of maximal TF induction to proceed. We increase 6.2-fold to reflect the altered maximal uptake from thrombin presence, hence = 6.2 = 26.23 1/min [46].

Synergistic effects of Thrombin/TNF
: To reflect the observed 6.2-fold reduction in TF induction due to TNF when thrombin is not present, we set 4 = 6.2ℎ , and = 1 [27].

Synergistic effects of TNF /IFN-:
We also capture the synergistic effects of TNF /IFN-on TF induction. IFN-has not been shown solely induce TF in vitro; however, IFN-modulates TF induction due to . This synergism of TNF /IFN-on TF induction is reflected in the first term of Equation (13) with Michael-Menten form.
Blood clotting severity parameter : We inform this parameter which appears in Equation (17) using in vitro experimental data [23]. To proceed, we first assume the following: 1) thrombin concentration is constant throughout the course of the experiment and 2) thrombin is spatially homogeneous in the artery. Writing Equation (17) Figure 4a in [23].

Initial conditions:
The initial conditions for model (1) are given in Supplementary Table S2. For the initial influenza virus dose and target epithelial cell count we use estimates and experimental data reported in [22]. Blood coagulation initial conditions were chosen to reflect in vitro experimental levels [12]. Realistically, circulating prothrombin and antithrombin levels vary person-to-person. High prothrombin levels have been correlated with a risk of arterial and venous thrombosis and are also affected by genetic conditions [20]. We capture these individual effects in our study by incorporating stochasticity in the initial values of prothrombin and thrombin for relevant simulations. In particular, we vary each initial condition by 20% of their baseline value. We begin simulations with model (1) at the disease free equilibrium, with an initial viral load, to reflect an infection of an individual at baseline status. We discuss these details in section 4.2-4.3 of the main text. The initial conditions are shown in Supplementary  Table S2.
Threshold level of TNF to induce TF production: Here we use clinical data to estimate the threshold . During a 7-day window within initial influenza infection, an elevated risk of acute myocardial infarction (AMI) has been observed [1]. Among 364 hospitalizations for AMI in 332 patients with laboratory-confirmed influenza, 20 occurred during 7-day post-infection risk window and the remaining 344 occurred outside of the risk window. The weekly admissions per week rates were 3.3 to 20 (incidence ratio 6.05) [1]. To estimate , we attribute these AMIs to blood clot formation; hence, given an influenza infection an individual has a probability 20/364 of experiencing an AMI within 7 days. To ensure the clotting frequency given by model (1) is consistent with this clinical data, we use a minimization process to estimate .
The process to find optimal was conducted as follows. First, let be the observed probability of AMI among study patients in [1] following infection, i.e., = 20 364 . We sample model parameters and initial conditions from their respective ranges (Supplementary  Table S3) to generate 364 parameter sets using Latin Hypercube Sampling. Now, let ( ) = 364 , where is the number of blood clot events in 364 model runs from all parameter sets, for a given threshold value of . We then search for an optimal such that we minimize the difference between model simulation and clinical observations. This process can be described mathematically as follows.
Let be the probability of AMI in 7 days. Find by solving the following minimization problem: We solve this numerically using Matlab's fminsearch function. The result is a threshold TNF level = 27.36 pg/ml. For more details about Latin Hypercube Sampling and optimization process, see [53].

S1.3 Rationale for IFN-I units
We have adopted the modelling construction of IFN-I kinetics from a prior work [22]. In this work, IFN-I was not directly observed and was left dimensionless. In the present study, we adopt this choice in units.