key: cord-0937361-d38q164g
authors: Fadai, Nabil T.; Sachak-Patwa, Rahil; Byrne, Helen M.; Maini, Philip K.; Bafadhel, Mona; Nicolau, Dan V.
title: Infection, inflammation and intervention: mechanistic modelling of epithelial cells in COVID-19
date: 2021-02-17
journal: Journal of the Royal Society, Interface
DOI: 10.1098/rsif.2020.0950
sha: 71367dad067fdcff96717ce8443423abf48f498b
doc_id: 937361
cord_uid: d38q164g

While the pathological mechanisms in COVID-19 illness are still poorly understood, it is increasingly clear that high levels of pro-inflammatory mediators play a major role in clinical deterioration in patients with severe disease. Current evidence points to a hyperinflammatory state as the driver of respiratory compromise in severe COVID-19 disease, with a clinical trajectory resembling acute respiratory distress syndrome, but how this ‘runaway train’ inflammatory response emerges and is maintained is not known. Here, we present the first mathematical model of lung hyperinflammation due to SARS-CoV-2 infection. This model is based on a network of purported mechanistic and physiological pathways linking together five distinct biochemical species involved in the inflammatory response. Simulations of our model give rise to distinct qualitative classes of COVID-19 patients: (i) individuals who naturally clear the virus, (ii) asymptomatic carriers and (iii–v) individuals who develop a case of mild, moderate, or severe illness. These findings, supported by a comprehensive sensitivity analysis, point to potential therapeutic interventions to prevent the emergence of hyperinflammation. Specifically, we suggest that early intervention with a locally acting anti-inflammatory agent (such as inhaled corticosteroids) may effectively blockade the pathological hyperinflammatory reaction as it emerges.

SARS-CoV-2, the causative viral agent for COVID-19 illness, has infected tens of millions of people and caused many deaths worldwide. As SARS-CoV-2 is a novel coronavirus, with ensuing COVID-19 a new illness, knowledge transfer has been piecemeal and only preliminary data exist for SARS-CoV-2 [1] . Treatment thus far has focused on severe COVID-19 [2] or the production of vaccines to prevent transmission [3, 4] . However, as it has been recognized that COVID-19 produces a 'runaway' train pattern of hyperinflammation [5, 6] , it can be postulated that early intervention with an anti-inflammatory, such as inhaled corticosteroids (ICS), could contain the illness [7] . This hypothesis is supported further by the clinical finding that chronic respiratory diseases (such as asthma and chronic obstructive pulmonary disease, COPD) are not the most common co-morbidities seen in COVID-19 [8] , indicating that ICS, a common treatment for asthma and COPD, could provide a protective response.

The majority of COVID-19 patients develop a respiratory illness [9] , with the primary site of infection being the airway epithelium where the ACE2 receptors, necessary for SARS-CoV-2 infection [10] , are in abundance. Furthermore, the emerging literature now contains evidence [7] that supports a conceptual model (figure 1) in which a viral infection triggers an immune response, which, in a minority of patients, is an immune-driven hyperinflammatory process [5] , leading to an acute respiratory distress-like syndrome, as well as causing direct damage by other (e.g. vascular) means [11] .

As we gain further understanding of SARS-CoV-2, it is clear that mechanistic investigations alone may not be sufficient and that mathematical models of a complex biological system are needed [12] . Specifically, mathematical models describing lung epithelial inflammation in COVID-19 may provide insight not accessible by experimental work alone [12] . Such approaches have previously been used to investigate the onset and resolution of inflammation [13, 14] , identify key parameters in the inflammatory process and better understand the interplay between inflammation and cell migration in epithelial tissue [15] . The urgent need better to understand COVID-19 pathology, combined with the necessarily partial and piecemeal experimental evidence emerging from clinics and laboratories, makes mathematical modelling a natural approach for hypothesis generation and testing, as well as designing treatment strategies.

We hypothesize that in at-risk patients, clinical deterioration is a combination of direct damage from the virus and a resulting uncontrolled pulmonary inflammatory process; we develop a mathematical model to explore these phenomena. In particular, we present a five-component mathematical model that describes inflammation kinetics, viral infection, and novel mechanisms that couple their dynamics. Using this model, we (i) investigate which qualitative differences in inflammation levels can be observed in patients infected by SARS-CoV-2, and (ii) determine how ICS treatment can be employed to reduce inflammation in the lung epithelium, thereby preventing the onset of hyperinflammation and other severe damage associated with high levels of pro-inflammatory mediators.

Using the current understanding of COVID-19 [5,6], we consider a mathematical modelling framework that describes the temporal evolution of five species at early timescales of SARS-CoV-2 viral infection: pro-inflammatory mediators (e.g. cytokines, interferons, C), recruited immune system cells (e.g. lymphocytes, T-cells and macrophages, M), freelymoving SARS-CoV-2 virus (V), cells susceptible to viral infection (e.g. epithelial lung cells, S), and infected cells (I). We refer to this five-species mathematical model as the MVSIC model. These species interact via a network of physiological responses relating to the localized infection and inflammation of epithelial lung cells (see figure 2 for a model schematic). Unlike previous models describing inflammation or viral infection (c.f. [13, [16] [17] [18] ), the MVSIC model incorporates the novel additional feature of coupling inflammation kinetics to viral infection.

We assume that each species will have a natural death/ clearance rate (process A in figure 2 ) and their dynamics are coupled via the following mechanisms. Firstly, we assume that susceptible epithelial lung cells follow logistic growth in the absence of infection (process B in figure 2) and become infected when free virus enters and reproduces inside it ( process D in figure 2 ). Secondly, we assume that infected cells produce several copies of free virus; these new viruses leave the infected cell without destroying it ( process E in figure 2 ). Thirdly, we assume that susceptible and infected cells can both 'signal' for cytokine recruitment [5, 6] after detecting free virus ( process F in figure 2 ); the infected cell signal may be smaller or even repress the susceptible cell's signals. This cell signalling mechanism is represented by rate-limited kinetics [19] , which assumes a sigmoidal Hill-like reaction term [20] , multiplied by the amount of signalling cells present (see §2.1 for further details). Fourthly, we assume that susceptible cells can also directly recruit (local) immune cells, bypassing the cytokine pathway, via ratelimited kinetics ( process G in figure 2 ). While this direct pathway allows faster recruitment of immune cells, there are fewer immune cells locally available. Therefore, we also assume that cytokines are recruited via immune system cells ( process H in figure 2 ) and vice versa ( process J in figure 2), creating a positive feedback loop. Finally, we assume that immune cells can consume infected epithelial cells via phagocytosis and remove SARS-CoV-2 virus ( process K in figure 2 ). Figure 1 . A simplified model of the potential inflammatory process in COVID-19 and potential therapeutic intervention. A rising viral load causes rising inflammation, which causes damage but is normally kept in check by a feedback mechanism. In COVID-19 and other respiratory infections, this mechanism appears to fail, leading, in a minority of patients, to a 'runaway train' pattern of hyperinflammation. Separately, viral entry into the vasculature leads to direct damage in the lung and in distal organs. Potential intervention early on with an anti-inflammatory agent targeted at lung epithelial cells may help to re-establish this check on hyperinflammation.

royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 18: 20200950

We are not only interested in the underlying physiological mechanisms driving inflammation in epithelial lung cells, but also the effect of various intervention strategies that can reduce viral infection and/or severe cases of inflammation. Therefore, we also consider two additional pathways associated with therapeutic intervention (green dashed lines in figure 2 ). One intervention pathway is the reduction of viral infection of susceptible cells (process L in figure 2 ), as might appear in retroviral therapies or vaccines; this therapy is represented by a reduction in the viral infection rate. The other intervention pathway that we will consider is the accelerated clearance of pro-inflammatory cytokines by anti-inflammatories ( process N in figure 2 ), such as ICS [21] .

We now propose a mathematical formulation of the network diagram shown in figure 2 associated with our five-species model, while incorporating the various physiological mechanisms outlined in the previous section. The temporal evolution of the five species outlined in the MVSIC model is described using the following system of ordinary differential equations (ODEs):

The aforementioned rate-limited mechanisms, associated with signalling and recruitment, are represented by the rational functions in V appearing in equations (2.4) and (2.5). All model variables and parameters appearing in this dimensional version of the MVSIC model are described in table 1. In particular, we note that our intervention parameters associated with reducing the viral infection rate and increasing the pro-inflammatory cytokine clearance rate are denoted as ϵ and ϕ, respectively. While there are a large number of (dimensional) parameters present in the MVSIC model, certain combinations of these parameters can exhibit the same qualitative features to one another if particular non-dimensional parameter groupings remain unchanged. As a result, we now examine a nondimensionalized version of the MVSIC model to better understand how the qualitative features of the model vary with respect to these non-dimensional parameter groupings. By rescaling the dimensional variables as

we obtain, after dropping hats, the non-dimensional system royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 18: 20200950

where

(2:12)

Since δ 1 , the growth rate of epithelial lung cells, typically evolves on the timescale of days [22] , which is a similar timescale to early stages of SARS-CoV-2 infection, this is an appropriate choice of timescale to non-dimensionalize the MVSIC model. Finally, the non-dimensional MVSIC model (2.7)-(2.11) also includes initial conditions for the five species to reflect that initially, there is a small amount of virus present, no susceptible cells have been infected yet, and no inflammatory mediators or immune cells are present at the site of infection:

For simplicity, and minimal loss of generality, we will assume from here onwards that both Hill power coefficients, n and m, are equal to 2. These power coefficients represent a slower uptake of recruitment at low concentrations, but a larger response at moderate and high concentrations [19, 20] . Furthermore, this choice of n, m > 1 agrees with the model description of having a faster local recruitment for small viral loads, but is surpassed at moderate viral loads by cytokine-related signalling pathways for immune cell recruitment.

One of the objectives for developing the MVSIC model is to create a mechanistic model capable of describing five main clinical features observed in patients infected with SARS-CoV-2: (i) healthy patients that clear the virus without experiencing symptoms, (ii) asymptomatic patients that do not clear the virus but do not experience symptoms, (iii-v) patients experiencing mild, moderate, or severe levels of inflammation, the latter of which are also incapable of clearing the virus on their own, resulting in a hyperinflammatory state. These traits are not only based on the percentage of susceptible cells present, but also on the levels of proinflammatory mediators in the system. In other words, the virus-free and mild inflammatory states may look similar in terms of percentage of healthy epithelial cells, but are distinguished by the differing levels of pro-inflammatory mediators present. As such, we will define three main categories of qualitative features of the MVSIC model, based on percentage of healthy susceptible cells: (a) virus-free/ mild inflammatory state, (b) asymptomatic/moderate inflammatory state, and (c) severe inflammatory state. Categories (a) and (b) are further subdivided between virus-free or mild inflammation (asymptomatic or moderate inflammation) based on the absolute magnitude of the cytokine levels. However, it is important to note that the dimensional quantity of pro-inflammatory mediators and immune cells present could indeed be large or small, depending on the parameters Mathematically, there are two key features that distinguish these three qualitative states: the absence or presence of an infectious steady state, and if this infectious steady state is stable or unstable. When the infectious steady state does not exist, the system settles to the healthy steady state, whereby the virus is cleared by the immune system and susceptible cells are replenished to their full quantities (S = 1). We will now examine additional stability features of these equilibria.

A question that naturally arises from the MVSIC model is what qualitative features do we expect to observe for a given set of parameter values. To answer this question, we examine the existence and stability of various steady states, or equilibria, present in the MVSIC model. A steady state is a combination of constant values (S*, I*, V*, M*, C*) for which the right hand sides of (2.7)-(2.11) are all zero. From (2.8) and (2.9), any equilibria of the MVSIC model must have that I* = V* = M* = C* = 0; we will refer to this equilibrium as the zero state, which can be shown to always be unstable. Alternatively, if S* = 1 − κ(1 − ϵ)V*, then (V*, M*, C*) satisfy the following nonlinear system of equations:

(2:17)

Notably, if V* = 0 in (2.15), this implies that M* = C* = I* = 0 and S* = 1; this equilibrium represents the healthy steady state. Otherwise, we obtain a new equilibrium in which V*, M*, C*, I*, S* are all strictly positive; we refer to this equilibrium as the infectious steady state. Furthermore, it can be shown (electronic supplementary material) that the infectious steady state can only exist when

In other words, when susceptible cells become infected faster than infected cells being cleared from the system and the amount of free virus is sustained above a threshold quantity, the system settles to an infectious steady state. It can be shown that this infectious steady state, when (2.18) holds, is unique (electronic supplementary material). Additionally, the healthy steady state is unstable only when (2.18) holds, as this inequality corresponds to when all the eigenvalues of the Jacobian of (2.7)-(2.11) at the healthy steady state no longer have negative real parts. Therefore, the system can only clear the virus and inflammation responses completely if and only if the infectious steady state does not exist. In addition to the existence of a unique infectious steady state, we can also numerically compute where this infectious steady state becomes unstable. As shown in figure 4 , the infectious steady state undergoes a Hopf bifurcation in parameter space, whereby two of the eigenvalues of the Jacobian of (2.7)-(2.11) at the infectious steady-state transition from having a negative real part to having a positive real part. This transition causes the infectious steady state to be locally unstable and, as seen in the severe inflammation state in figure 3 , exhibit large oscillations in all species concentrations.

We also note that the ICS intervention parameter associated with pro-inflammatory cytokine clearance, ϕ, does not appear in (2.18) . While increasing cytokine clearance cannot remove the infectious steady state, we will show in the next section that this parameter is linked to a reduction in the pro-inflammatory mediators, thereby eliminating a potential hyperinflammatory response.

We now examine how two key pathways in the MVSIC model, related to the viral infection rate and the pro-inflammatory cytokine clearance rate, are altered in the presence of a medical intervention, such as ICS. One intervention strategy is to reduce the rate at which susceptible cells become infected, as might be achieved by a vaccine. In the MVSIC model, this therapy corresponds to the parameter ϵ, which does indeed drive severe cases of inflammation down to the virus-free state as ϵ increases, as shown via the inequality (2.18) . However, in the absence of a vaccine option for SARS-CoV-2, we must also consider other treatment strategies that reduce inflammation, even if it does not completely eliminate the virus from a patient.

Another intervention strategy is to increase the clearance rate of pro-inflammatory cytokines, thereby reducing the amount of inflammation present near infected cells. Since inhaled corticosteroids (ICS) act directly on epithelial lung cells, this therapy corresponds to increasing the parameter ϕ without any time delay. While ICS intervention cannot remove the infectious steady state, it can reduce large fluctuations of pro-inflammatory cytokines and the resulting damage that occurs in cases of severe inflammation. In figure 5a, we observe that the pro-inflammatory cytokine levels in the severe inflammation parameter regime (table 2) are reduced to about 25% when ϕ is increased from 0 to 3. Therefore, even in the case where ICS do not decrease the viral infection rate, they are still able to significantly reduce inflammation and augment the natural clearance rate of pro-inflammatory mediators. Indeed, by varying ϕ in magnitude, we see in figure 5b that the height of the proinflammatory cytokine spike decreases as ϕ increases. Relative to no intervention, i.e. ϕ = 0, we note that the spike's height decreases approximately by the relative factor F(ϕ) = 1/(1 + ϕ), which agrees with the numerical solution of the spike decrease factor with high accuracy (figure 5b). In other words, if ICS augment the natural cytokine clearance rate by a multiplicative factor of 1 + ϕ, then the height of cytokine spikes are reduced by approximately 1/(1 + ϕ) (e.g. doubling the effective clearance rate halves the spike height).

We have developed the first mathematical model describing the dynamics of inflammation arising in epithelial lung cells infected with SARS-CoV-2. This model, named the MVSIC model, incorporates a network of mechanistic and The pathological mechanisms in COVID-19 illness are still being elucidated and very much an active topic of investigation. However, it is recognized that patients infected with SARS-CoV-2 have high levels of pro-inflammatory mediators [25] , especially those that have severe illness. In the lung itself, pro-inflammatory monocyte-derived macrophages are abundant in the bronchoalveolar lavage fluid from patients with severe COVID-19 [26] . These, and other, findings point to a hyperinflammatory state in severely ill patients, which is believed to be linked to poor or fatal outcome, with a clinical trajectory that resembles acute respiratory distress syndrome (ARDS). Therefore, understanding, and preventing, dangerously high levels of inflammatory mediators present in patients with COVID-19 would appear to be crucial. Our model is consistent with these immunological findings; furthermore, this analysis points to potential therapeutic interventions to prevent the emergence of hyperinflammation (e.g. UK clinical trial NCT04416399 [7] ). Specifically, we suggest that an early intervention with a locally acting (i.e. targeted at inflamed epithelial lung cells) anti-inflammatory agent may effectively lead to blockade in the 'runaway train' inflammatory reaction. Inhaled corticosteroids or other cytokine-reducing medications, as recently suggested [7] , are candidates for such agents.

Our model is based on a simplified picture of the inflammatory pathways involved and therefore has limitations. It does not account comprehensively for every mechanism and cell class involved in SARS-CoV-2 infections and damage. We do not, for instance, distinguish between T-cells, macrophages, lymphocytes, and other immune cell types, instead focussing on the inflammatory mechanism as a whole. Additionally, in severe cases of COVID-19, recent evidence indicates that the virus permeates into the endothelium tissue of the lungs, circulating through the blood and causing distal organ damage by vascular and other means [27] , while here we do not consider this direct/ non-inflammatory avenue of viral damage. The timescale associated with SARS-CoV-2 permeating into the endothelium is far longer than local epithelial infection kinetics. Thus, the MVSIC model dynamics are relevant chiefly during the earlier stages of COVID-19 illness, when inflammatory processes dominate the clinical picture.

Finally, while certain immune responses, such as phagocytosis, produce inflammatory mediators, other immune responses, such as T-cells removing virus, do not. Thus, more detailed mathematical models are likely to be needed to account for the complex interplay between differing classes of immune cells; however, the MVSIC model is a first attempt to understand the biological complexity, which will require prospective validation. Enhancing our knowledge using these mathematical models is expected to be valuable in designing more sophisticated and patient-specific intervention strategies.

Data accessibility. All Matlab codes used to generate figures 3-5 are available at https://github.com/nfadai/MVSIC2020 helpful comments. Figure 5 . Intervening inflammation using inhaled corticosteroids. The severe inflammation parameter regime in table 2 is used prior to intervention (0 < t < 20). (a) At t = 20 (black arrow), ϕ is increased from 0 to 3 (red solid curve), reducing pro-inflammatory cytokine levels. The C(t) trajectory without intervention (ϕ = 0) is shown as a dark red dashed curve. (b) The height of the pro-inflammatory cytokine spike, relative to no intervention, decreases as ϕ increases. The numerically computed spike height decrease is shown in solid black, while the approximate decrease factor F(ϕ) = 1/(1 + ϕ) is shown in dashed blue.

SARS-beginning to understand a new virus

Dexamethasone in hospitalized patients with COVID-19-preliminary report

An mRNA vaccine against SARS-CoV-2-preliminary report

Safety and immunogenicity of the ChAdOx1 nCoV-19 vaccine against SARS-CoV-2: a preliminary report of a phase 1/2, single-blind, randomized controlled trial

2020 Pathological inflammation in patients with COVID-19: a key role for monocytes and macrophages

The trinity of COVID-19: immunity, inflammation and intervention

Inhaled corticosteroids in virus pandemics: a treatment for COVID-19?

Prevalence of comorbidities in the novel Wuhan coronavirus (COVID-19) infection: a systematic review and meta-analysis

Clinical course and risk factors for mortality of adult inpatients with COVID-9 in Wuhan, China: a retrospective cohort study

Structure of the SARS-CoV-2spike receptor-binding domain bound to the ACE2 receptor

Extrapulmonary manifestations of COVID-19

Mathematical modeling of complex biological systems: from parts lists to understanding systems behavior

The resolution of inflammation: a mathematical model of neutrophil and macrophage interactions

Identifying important parameters in the inflammatory process with a mathematical model of immune cell influx and macrophage polarization

Chronic TNFα-driven injury delays cell migration to villi in the intestinal epithelium

Mathematical modelling of atherosclerosis as an inflammatory disease

Modelling viral and immune system dynamics

Human immunodeficiency virus infection: from biological observations to mechanistic mathematical modelling

First-order stochastic cellular automata simulations of the Lindemann mechanism

The Hill equation and the origin of quantitative pharmacology

The cellular kinetics of lung alveolar epithelial cells and its relationship with lung tissue repair after acute lung injury

A quantitative model-independent method for global sensitivity analysis of model output

Variance based sensitivity analysis of model output. Design and estimator for the total sensitivity index

A dynamic COVID-19 immune signature includes associations with poor prognosis

Single-cell landscape of bronchoalveolar immune cells in patients with COVID-19

2020 Coronaviruses pathogenesis, comorbidities and multi-organ damage-a review

Acknowledgements. We thank the three anonymous referees for their