© The authors. Published by EDP Sciences, 2020

1 Introduction

On March 22, 2020 according to the official data presented at the Worldometer [19] for COVID-19, the number of infected individuals exceeded 10 thousand in seven countries. The respective growth rate of the number of infected is shown in Table 1.

The second column represents the number of days needed to pass from 100 to 1000 infected, the third column from 1000 to 10000. If the number of days is the same in the two columns, then the growth is exponential, as it is predicted by conventional epidemiological models (SIR) for the beginning of epidemic when the number of uninfected (susceptible) individuals can be considered as approximately constant. In Italy, Iran, and France, the number of days in the third column is larger than in the second column. In Germany this difference is small at the moment, so it is early to make any conclusions. In USA and Spain the number of days is the same and the growth rate is slightly larger in Spain. A similar scenario was observed later in other countries.

These data should be taken with caution because they can be strongly influenced by the number of effectuated tests. According to some media reports, there is a shortage of testing facility and/or shortage of testing kits in many countries, and the accurate number of infected individuals can be much larger than the number reported so far [7].

Dynamics of the number of infected individuals can be influenced by various factors, in particular, by their spatial distribution and nature of social mixing. Appearance of new focuses of infection can essentially change the growth curves. For example, at the whole world level, development of the epidemic in Europe started when it was mostly at decaying mode in Asia. Instead of the standard growth-decay curves, this leads to growth-decay-growth curves. Some elements of this more complex dynamics where the number of newly reported daily cases decreased during several days was also observed in some countries.

If we suppose that the distribution of infected individuals is uniform over space, then the SIR model predicts the exponential growth at the beginning of epidemic. Some time later, the growth rate decreases due to the decrease in the number of susceptible individuals, and the number of infected individuals decays at the end of epidemic. There are recent works with various modifications of the epidemiological models trying to take into account some specific features of the coronavirus pandemic [6, 9, 18].

Analysis of the data on seasonal influenza shows that dynamics of epidemic spread can include the period of accelerated growth (Fig. 1). In the beginning of epidemic, the growth rate is exponential. At the next stage it continues to grow exponentially but with an increased growth rate. It can be considered as two different exponential functions where the second one has a larger exponent than the previous one. A possible explanation of this effect is related to the presence of sub-populations with different response on the infection. Simplifying the situation, we can consider two sub-populations, one of them with a weak immune response and another one with a strong immune response. The strength of immune response is determined by the rate of production of antigen specific immune cells. We will define it below mathematically. At the beginning of epidemic, infection propagates among the first sub-population. When the epidemic reaches certain level, the second sub-population becomes also involved. This increase of the number of susceptible individuals result in the acceleration of the growth rate.

Transition between the two stages of epidemic growth can be related to the number of infected individuals. If it is small enough, then a healthy individual rarely meets infected individuals, and the initial viral load remain small. People with a strong immune system can eliminate infection without developing disease symptoms, while people with a weak immune system can fall sick. At the later stages of epidemic, when the number of infected individuals is significantly large, an uninfected individual can cross several infected individuals during a short period of time. Therefore, the initial viral load becomes larger, and people with a strong immune response can also become ill.

In order to describe this two stage epidemic growth, we develop in this work an immuno-epidemiological model combining the immune system at individual level and the epidemic at population level. This model is generic but it can take into account specific features of a particular infection or the details of epidemiological situation (like quarantine, social distancing, self-isolation, etc.). In the next section, we will discuss some epidemiological data on influenza epidemic. Section 3 is devoted to a model of immune response where we show how people with weak and strong immune system can respond on a given initial viral load. Immuno-epidemiological model is studied in Section 4.

Figure 1 Data on seasonal influenza in Europe in 2019–2020 [5]. The left graph shows the number of influenza-positive specimens from non-sentinel sources (blue line) and its approximation by an exponential (red line). The right graph shows the same curves in the logarithmic scale.

2 Two stage epidemic

There are large amount of data available for the seasonal influenza epidemic. A typical graph of the number of cases is shown in Figure 1 (left, blue curve) [5]. Growth rate of the number of positive cases between weeks 1 and 6 is exponential. The red curve shows the approximation of the blue curve with the function 20.5 exp(0.2x), where x is the number of weeks. The same curves are shown in Figure 1 (right) in the logarithmic scale. The red curve here is a straight line. The weeks 7–10 show different growth rate with some acceleration in growth. It also correspondsto an exponential growth but with a different exponent. Some additional data are presented in Appendix A.

Conventional epidemiological models like SIR (Susceptible, Infected, Recovered, see, e.g., [12]) and their numerous variants describe the density of infected individuals I by the typical equation (2.1)

where S and I are the densities of susceptible and infectedindividuals, and λ, μ are positive constants. The first term in the right-hand side of this equation characterizes the appearance of new infected individuals due to their contact with susceptible. The second term corresponds to the decrease of I due to the recovery or death of infected individuals.

In the beginning of the infection spreading, the number of infected and recovered individuals are much less than the number of susceptible, so that we can approximate S by a constant S ≈ S₀. Using this approximation, we obtain a linear differential equation with constant coefficients. Its solution is given by the function I(t) = I₀ exp(αt), where I₀ is the number of infected individuals at the initial moment of time, and α = λS₀ − ν. If α > 0, then I(t) grows exponentially. The same condition can be written as , where is the basic reproduction number.

Thus, according to equation (2.1), growth of the number of infected individuals is exponential in the beginning of epidemic, and it slows down later when S decreases. Exponential growth in the beginning corresponds to the data in Figure 1. However, it does not describe accelerated growth during the second stage of epidemic. A possible explanation of this acceleration is that in the beginning the epidemic, disease spreads mainly among the people with weak immune response. At the second stage, people with strong immune response become also exposed. In order to investigate this hypothesis, we consider a model of immune response in the next section.

3 Model of immune response

We considera qualitative model of immune response suggested in [2]: (3.1) (3.2)

where v and c are the concentrations of viruses and immune cells within the organism. The first term in the right-hand side of equation (3.1) describes the virus multiplication rate. This logistic term is proportional to the virus density v and to the normalized density of uninfected cells (1 − v). The second term characterizes virus death due to the immune response. This term is proportional to the virus density and to the density of immune cells c. The first term in the right-hand side of equation (3.2) describes multiplication of immune cells and the second term their mortality. In general, both of them depend on virus concentration. Clonal expansion of immune cells is stimulated by antigen (virus) if its concentration is not very large. For large virus concentration, the proliferation rate of immune cells can decrease. Therefore, the function p(v) is non-negative and increasing with saturation or increasing for small v and decreasing for large v. The death rate of immune cells can also depend on the concentration of viruses. In this case it is a positive increasing function. It can be also considered as constant if this effect is independent of virus density.

The nullclines of system (3.1), (3.2) are given by the expressions: v = 0 and c = k(1 − v)∕σ; c = 0 and c = f(v), where f(v) = 1 − h(v)∕p(v). Figure 2 qualitatively shows possible phase portraits of this system. There are six stationary points in the left phase portrait and four in the right one. If there is a small initial virus load v₀ and a smallnumber of antigen specific immune cells (new infection, no memory cells), then the trajectory of this system will start from vicinity of the unstable stationary point (0, 0) and move towards one of the other stationary points. Summarizing the behavior of solutions, we observe that if the initial viral load v₀ is sufficiently small, then the trajectory approaches the stationary point with complete curing (Fig. 2, right) or with a weak persistent infection (Fig. 2, left). If the initial viral load is large enough, then the solution converges to a stationary point with a large virus concentration and low concentration of immune cells. It can characterizea strong chronic infection or a lethal outcome. We present some results of numerical simulations in Appendix B.

In order to get an explicit solution, we simplify system (3.1), (3.2) with the following assumptions:

• (A1)

  functions p(v) and h(v) are constant. This means that the clonal reproduction of immune cells fully responds to the antigen as soon as there is a small amount of virus in the organism, and mortality rate of the immune cells does not depend on virus concentration;

• (A2)

  there is no limitation on the production of viruses and immune cells by available resources, i.e., the logistic terms in equations (3.1), (3.2) are replaced by the linear terms.

Under these assumptions, instead of system (3.1), (3.2) we consider the following system of equations: (3.3) (3.4)

We find from equation (3.4), c = c₀e^μt, where c₀ is the number of initial antigen specific immune cells, and μ = p − h. Assuming that μ > 0, we obtain exponential growth of immune cells. Substituting the solution for c into equation (3.3), we find the solution for v as

where v₀ is the initial viral load. The derivative of this function

is positive at t = 0 since the initial number c₀ of antigen specific cells is very small. Therefore, v(t) grows at the beginning of infection development. This derivative equals 0 at time when the maximal virus concentration is reached. After that, it decreases and converges to 0. The maximal virus concentration is as follows:

where w = k∕(c₀σ), s = (p − h)∕(c₀σ). Parameter w corresponds to the dimensionless virus multiplication rate. It characterizes the virulence of infection. Parameter s characterizes the strength of the immune response.

The maximal level of viral concentration is determined by the initial viral load v₀, the virulence of infection w, and the strength of immune response s. For a given virulence of infection w > 1, the maximal viral load increases for large initial viral loads and decreases for strong immune response. If we set W = wlnw − w + 1, then v^* = v₀ exp(W∕s).

Suppose that there are two levels of immune response s₁ and s₂, with s₁ < s₂. The first one corresponds to the weak immune response and s₂ to the strong immune response. For a given initial viral load v₀, the corresponding maximal virus concentrations are and , where . If virus concentration exceeds some critical value, which is determined by a proportion of infected tissue, then disease symptoms are developed, and the person is considered as infected and fall sick. If this critical level of virus concentration v_c is such that , then people with weak immune response become ill while people with strong immune response remain healthy. The same property is preserved in the full model, where s₁ < s₂ leads to for the same v₀ (see Appendix B).

However, if the initial virus load increases with the progression of epidemic, then also increases and it can overpass the critical level. We consider this question in the next section.

Figure 2 Two examples of phase portrait of system (3.1), (3.2) depending on the function f(v) = 1− h(v)∕p(v). It can have six (left) or four (right) stationary points. The trajectory with the initial condition (v0, c0 ) and sufficiently small viral load v0 can go either to the stationary point with complete virus elimination (right) or to a weak persistent infection (left). If the initial viral load is sufficiently large, then the trajectory approaches the stationary point with large v and small c. It corresponds to a strong chronic infection or to death. Reprinted from [2].

4 Immuno-epidemiological model

We can now combine the epidemiological model and the immune response model making a link between the number of infected individuals I, the initial viral load v₀, and the number of susceptible individuals S. We suppose that v₀ = v₀(I) and S ≈ S(v₀) = S(v₀(I)). Hence, the number of susceptible becomes a function of the number of infected.

4.1 SIR model

Under the assumption that the number of susceptible individuals depends on the number of infected through the initial viral load, we obtain a closed equation for the number of infected: (4.1)

As before, let us assume that the total population consists of two sub-populations, one of them with weak immune response and another one with strong immune response. If the initial viral load v₀ is less than some critical value , then the maximal viral concentration in the case of strong immune response remains less than the critical level v_c which determines the appearance of the disease. This critical viral level can be determined from the equality . Hence

where S₁ correspondsto the sub-population with weak immune response and S₂ to the total population. Under this assumption, the explicit dependence v₀ (I) is not needed. We suppose that the viral load is an increasing function such that for I < I_c, and for I ≥ I_c. Then instead of equation (4.1) we can write: (4.2)

We can now find the solution:

Here α₁ = λS₁ − ν, α₂ = λS₂ − ν. Hence, the number of infected individuals grows with exponent α₁ till some critical time t_c, and with a larger exponent α₂ for t > t_c.

4.2 Quarantine model

Another epidemiological model was introduced in [18], in order to take into account removal of infected individuals to quarantine when the incubation period is finished and they manifest disease symptoms: (4.3)

Here I(t) is the density of latently infected individuals (during the incubation period), τ is the length of the incubation period. In this case we have three time intervals: 0 < t < t_c where t_c is such that I(t_c) = I_c; t_c < t < t_c + τ, and t > t_c+ τ. We assume S(t) = S₁, S(t −τ) = S₁ for 0 < t < t_c, and the solutionof equation (4.3) is as follows:

where I₀ is the initial number of latently infected individuals, and μ₁ satisfies the equation μ₁ = kS₁(1 − e^−μ₁τ). This equation has a positive solution if the basic reproduction number is greater than 1. Let us recall that characterizes the number of newly infected individuals during the incubation period. The value t_c can be found from the equality I(t_c) = I_c, .

During the second time interval t_c < t < t_c + τ, we have S(t) = S₂ and S(t − τ) = S₁. The solution of equation (4.3) writes:

where μ₂ satisfies the equality μ₂ = λS₂(1 − θe^−μ₂τ), θ = S₁∕S₂. Since θ < 1, this equation has a positive solution andhence μ₂ > 0.

If t > t_c + τ, then S(t) = S₂, S(t − τ) = S₂,

where I₁ = I_ce^μ₂τ, and μ₃ satisfies the equation μ₃ = kS₂(1 − e^−μ₃τ). This equation has a positive solution if the basic reproduction number is greater than 1.

Let us note that μ₁ < μ₃ < μ₂. Therefore, the final growth rate is larger than the initial one but less than the intermediate.

5 Discussion

Numerous data on seasonal influenza epidemic shows that growth of the number of infected individuals can occur in two stages. Both of them are characterized by an exponential growth rate but the second exponent is larger than the first one. Growth accelerates when the number of infected reaches certain threshold. It should be noted that usually several virus strains are observed during influenza epidemic. All of them show similar growth pattern.

Such dynamics can not be captured by conventional epidemiological models like SIR (susceptible-infected-recovered). These models predict exponential growth at the beginning of epidemic when the number of susceptible S can be considered as approximately constant or too large compared to initial number of infected. During the disease progression, S decreases, as a consequence the growth rate of disease progression also decreases. The number of infected passes through a maximum and then decays to 0.

The goal of this work is to develop a minimal mathematical model which can describe accelerated growth of disease spread. Our main hypothesis is that this acceleration can occur if we take into account the presence of different sub-populations, namely, with weak and strong immune response. In the beginning of epidemic, only (or mainly) people with weak immune response are concerned while at the second stage people with strong immune response also become susceptible. This transition takes place when the number of infected individuals becomes significantly large, and an uninfected individual has enhanced probability of interactions with infected individuals. This increases the viral load received during a short period of time and leads to the disease propagation in the sub-population with strong immune response.

In order to analyze how the initial viral load influences the disease initiation and progression, we consider an immuno-epidemiological model describing the concentration of viruses and immune cells in the organism. The strength of immune response in the model is determined by the rate of multiplication of immune cells. We show that for the same initial viral load, an individual with a weak immune system falls ill while an individual with a strong immune system does not develop disease symptoms. However, if the initial viral load increases, then the individuals with strong immune response also start developing the disease.

The influence of the initial viral load on the disease progression and outcome was studied for various influenza epidemics and animal models [11]. It is shown that it increases the disease severity and related mortality. A dose-dependent response determined by the number of simultaneous contacts between susceptible and infectious individuals could determine the dynamics of the influenza pandemic in 1918 [13].

The combination of immunological and epidemiological models allows the description of the interplay between the population level and individual level [3, 4, 10, 14, 17]. We use here this approach to study a possible acceleration of epidemic growth rate.

We describe a generic situation applicable to different types of epidemics. During the influenza epidemic, the critical level of incidences when the growth accelerates is from 50 to 150 per week over 100,000 of the population (Figs. 1, A.1, A.2). Taking into account that the disease duration is about one week, then we get an estimate of the number of simultaneously infected individuals. The situation is different for the coronavirus infection because the disease duration is longer (about 4 weeks) and because the patients with manifested symptoms are put in quarantine. Therefore, we need to take into account here mainly the latently infected individuals during the incubation period. For a country with the total population 10 million and duration of the incubation period one week, the critical level of epidemic is about 1400 new daily cases. This estimate does not take into account the restrictions on the population movement (it does not affect the influenza epidemic).

The immuno-epidemiological model developed in this work is not pathogen specific. We suppose that it can be also applied to coronavirus infection since the strength of adaptive immune response plays important role in the disease progression [8]. It should be noted that in some cases excessive immune response can damage lung tissue and provoke important clinical consequences.

5.1 Limitations of the model

The main hypothesis of this work is that epidemic progression begins with a sub-population characterized by weak immune response and, at a later stage, it continues in the whole population. Though this assumption seems plausible, we do not have direct confirmation that this heterogeneity of the population plays a significant role during epidemic spread.

A particular choice of the epidemiological model is not essential here since all of them give exponential growth for the number of infected individuals. The exponential growth rate is proportional to the density of susceptible individuals. Once we suppose that effective density of susceptible can increase during epidemic because different sub-populations are involved one after another, we can obtain the increase of the growth rate.

The immunological model used in this study is very simplified. It does not take into account for the involvement of many different cell types, intracellular regulation, time delay in virus production and clonal expansion of immune cells, and some other relevant aspects. However, it captures the main features of the interaction of viral infection with the immune response. Its simplicity allows to couple it with the epidemiological model.

Let us note that system (3.1), (3.2) can also be explicitly solved if the functions p(v) and h(v) are linear (and the second assumption is satisfied).

Acknowledgements

The last author acknowledges the IHES visiting program during which this work was done. The work was supported by the Ministry of Science and Education of Russian Federation, project number FSSF-2020-0018, and by the French-Russian program PRC2307.

Appendix A Influenza epidemic data

We present here some typical data on influenza epidemic (Figs. A.1, A.2). The two stage growth can be observed.

Figure A.1 Data on seasonal influenza in UK [16]. The left graph shows the number of influenza-positive specimens from non-sentinel sources (blue line) and its approximation by an exponential (red line). The right graph shows the same curves in the logarithmic scale.

Figure A.2 Data on seasonal influenza in Belgium [15]. The left graph shows the number of influenza-positive specimens from non-sentinel sources (blue line) and its approximation by an exponential (red line). The right graph shows the same curves in the logarithmic scale.

Appendix B Solution of the immunological system

We present here some results of numerical simulations of system (3.1), (3.2) with the functions p(v) = (p₀ + p₁v)e^−p₃v, h(v)= h₁v, where p₀ =0.2, p₁ =5 (Fig. B.1) and p₁ = 10 (Fig. B.2), p₃ = 4, h₁ = 0.2. Other two parameter values are k = 1 and σ = 1.

Figure B.1 Trajectories on the phase plane of system (3.1), (3.2) in the case of weak immune response (left) and two sample plots for the growth of v(t) with time (right).

Figure B.2 Trajectories on the phase plane of system (3.1), (3.2) in the case of strong immune response (left) and two sample plots for the growth of v(t) with time (right).

The strength of immune response is characterized by the parameter p₁ which determines the rate of clonal expansion of immune cells in response to the antigen. A lesser value of p₁ correspondsto weaker immune response (Fig. B.1), a larger value of p₁ to stronger immune response (Fig. B.2).

There are 7 trajectories in each figure starting from two nearby initial points for some small values of v = v₀ (initial viral load) and c = c₀ (initial concentration of antigen specific cells). If the initial viral load is sufficiently small, then in both cases the trajectories converge to the virus-free stationary point corresponding to the curing. However, in the case of strong immune response the maximal and the total virus concentrations are less. Therefore, strong immune response can prevent the development of the disease while weak immune response may not be sufficient for this. If the initial viral load is sufficiently large, then weak immune response can lead to a chronic infection or to death while strong immune response leads to virus eradication.

References

All Tables

All Figures

	Figure 1 Data on seasonal influenza in Europe in 2019–2020 [5]. The left graph shows the number of influenza-positive specimens from non-sentinel sources (blue line) and its approximation by an exponential (red line). The right graph shows the same curves in the logarithmic scale.
In the text

Figure 2 Two examples of phase portrait of system (3.1), (3.2) depending on the function f(v) = 1− h(v)∕p(v). It can have six (left) or four (right) stationary points. The trajectory with the initial condition (v0, c0 ) and sufficiently small viral load v0 can go either to the stationary point with complete virus elimination (right) or to a weak persistent infection (left). If the initial viral load is sufficiently large, then the trajectory approaches the stationary point with large v and small c. It corresponds to a strong chronic infection or to death. Reprinted from [2].

In the text

	Figure A.1 Data on seasonal influenza in UK [16]. The left graph shows the number of influenza-positive specimens from non-sentinel sources (blue line) and its approximation by an exponential (red line). The right graph shows the same curves in the logarithmic scale.
In the text

	Figure A.2 Data on seasonal influenza in Belgium [15]. The left graph shows the number of influenza-positive specimens from non-sentinel sources (blue line) and its approximation by an exponential (red line). The right graph shows the same curves in the logarithmic scale.
In the text

	Figure B.1 Trajectories on the phase plane of system (3.1), (3.2) in the case of weak immune response (left) and two sample plots for the growth of v(t) with time (right).
In the text

	Figure B.2 Trajectories on the phase plane of system (3.1), (3.2) in the case of strong immune response (left) and two sample plots for the growth of v(t) with time (right).
In the text