Mathematical Modeling for CoronaVirus I: SI, SIR and SEIR Models

Outline

• ODE models
• Complex network models
• Statiscal models

In ODE models, divide the total population into several compartments and find ODEs between them.

SI model

Two compartments $S(t)$ and $I(t)$. Both of them are functions of time $t$.

• $S(t)$ Susceptible population;
• $I(t)$ Infected population.

Model

Conservation of population $S(t) + I(t) = N$, where $N$ denotes the total population.

with initial condition $S(0) = S_0, I(0) = I_0,$ and $I_0 + S_0 = N$.

Here $\beta$ is the rate of infection and $I/N$ is the fraction of infected people.

A compact representation is:

Explanation: $\frac{I}{N} S$ is the number of susceptible person will encounter an infected one. Multiply the rate of infection $\beta$ is the rate of change of infected person. The minus sign in the first equation is from the converstation of total population.

Simply add these two equations, we have $\frac{{\rm d}}{{\rm d} t}(S(t) + I(t)) = 0$ and thus $S(t) + I(t) = N$ for all $t$. Again this is the conservation of population. No death is included in the model.

Analytic solution

Substitute $S=N-I$ into the second equation to get a first order ODE

Its solution is

which is a standard logistic function and also known as the sigmoid function used in machine learning.

Solution to the ODE

$$x' - a x + bx^2 = 0, \quad x(0) = x_0.$$

Let $y = 1/x$ and obtain a linear ODE $y' + ay =b$ for which the exact solution is $y(t) = c e^{-at} + b/a$ where the constant $c$ is determined by the initial condition $y(0) = 1/x_0 = c + b/a$. Therefore the solution of (5) is

$$x(t) = \frac{1}{c e^{-at}+b/a}, \quad c = \frac{1}{x_0} - \frac{b}{a}.$$

The asymptotic of the solution is

A sad conclusion: eventually all people will get infected.

We plot several curves of $I(t)$ for different $\beta$ with other parameters fixed $N = 1000, I_0 = 1$. Although $\lim _{t \rightarrow+\infty} I(t)=N$ for all $\beta$, the derivative for different $\beta$ will have different peak at different time.

Numerical solution

For SI model, we can use the analytic formula but in general it is not easy to solve ODE systems analytically. The simplest numerical methods for solving ODE is the forward Euler method. Recall the definition of derivative

In forward Euler method, we chose a fixed $\Delta t$ and rewrite the ODE systems as the difference equation

To understand the meaning of (7), let us chose $\Delta t = 1$. Depending on the unit of the time, it means $1$ day, $1$ hour, or $1$ month. To fix the idea, we chose the unit as "day".

The meaning of equation for $I$ can be understood as: the number of infected people in the next day $t+1$ is equal today's $I(t)$ plus the new infected which is the susceptible $S(t)$ multiply the possibility of encounter an infected person the ratio of infected people $I/N$ and then multiply the infection rate $\beta$ .

Start from $I(0) = I_0, S(0) = N - I_0$, we can compute $I(1), S(1)$ and repeat the procedure to compute $I(t), S(t)$ for all later time.

The main drawback of the forward Euler method is the stability. If $\Delta t$ is too large, then the numerical solution could below up and not a faithful approximation of the true one.

Read: Numerical methods for ordinary differential equations for a short introduction and take Math 107 or Math 225.

Discussion

Due to the symmetry of the solution, we know when $I(t^*) = N/2$, the second derivative $I''(t^*) = 0$ . Solving this equation to get $t^* = \beta^{-1}\ln (N/I_0 -1)\approx \beta^{-1}\ln N$.

The first derivative of $I(t)$ is always positive which means $I(t)$ is a monotone increasing function but the speed of increasing acheives the peak at $t^*$.

This is almost true in related but more sophiscated models. So a rule of thumb for the predicition of the total infected number: when you observe the turning point (daily increasing number is decreasing), double the current number.

The fact all peopole will get infected is quite frustrating. What is the limitation of SI model? At least we didn't consider:

• Recovery. Infected person can get recovered from the illness.
• Death. Infected person can die.

To include the death, we introduce another parameter $\gamma$: the death rate. Then simply add a damping term to the second equation

Now adding these two equations will give a damping of the total population $\frac{{\rm d}}{{\rm d} t}(S + I) = -\gamma I \leq 0$.

We shall discuss SIR model to include the recovery component in the next section and continue to conisder the situation where the disease is not recoveryable. One such example is AIDS. Then SI model implies if we have some people get HIV positive, then eventually we all get it.

Obviously, this is not true. What's wrong here?

In ODE models, we assume people are full connected but this assumption is wrong for AIDS. To include the connectedness into the model, we should introduce the concept of complex network and use graph theory which will be discussed later on.

Exercise Try to analyze $\eqref{SIdeath}$ and discuss the behavior of the model for different parameters $\beta$ and $\gamma$ .

SIR model

We introduce one more compartments and consider SIR model.

• $S(t)$ susceptible population;
• $I(t)$ Infected population;
• $R(t)$ recovered population

Model

We still use parameter $\gamma$ to denote the rate of transition from infectious to recovered. Then the model becomes

with initial condition $S(0) = S_0, I(0) = I_0, R(0) = 0$, and $I_0 + S_0 = N$. Note that the first two equations are the same as $\eqref{SIdeath}$, where we interpret $\gamma$ as the death rate.

And the compact representation is

Again conservation of population can be obtained by summing all equations together. That is $S(t) + I(t) + R(t) = N$ for all $t\geq 0$.

Numerical solution

It is not easy to solve this system to get the analytic solution. But numerical simulation is still the same. We write out the forward Euler's method with $\Delta t = 1$.

which gives a better explanation of the model. For example, the last one means tomorrow's recovered people is equal to today's add the recovered ratio $\gamma$ multiply the number of infected people $I$.

In MATLAB, one can simply use ode45.

​x
function SIR(N,I0,beta,gamma,tend)
tspan = [0,tend];
S0 = N - I0;
R0 = 0;
y0 = [S0; I0; R0];
[t,y] = ode45(@(t,y) SIRfunc(t,y,beta,gamma,N), tspan, y0);
plot(t,y(:,1),'-',t,y(:,2),'m-',t,y(:,3),'g-','LineWidth',2);
legend('Susceptible','Infectious','Recovered')
title('SIR model')
end
​
function dydt = SIRfunc(~,y,beta,gamma,N)
    dydt = [-beta/N*y(2)*y(1); ...
             beta/N*y(2)*y(1) - gamma*y(2); ...
             gamma*y(2)];
    
end

Run SIR(1000,10,0.1,0.05,300) we get the following figure

Parameters for the figure: $N = 1000, I_0 = 10, \beta = 0.1, \gamma = 0.05$.

Observation

From the figure, we can conclude:

1. It will be under-control. The infected number will decrease to zero eventually.
2. A fixed fraction of the population becomes infected and eventually recovered and the rest population remains healthy.
3. There is a peak of the infected number.

Exercise

Prove there exists $S_{\infty}, R_{\infty}\in (0, N)$ so that $\lim _{t \rightarrow \infty} S(t)=S_{\infty}$ and $\lim _{t \rightarrow \infty} R(t)=R_{\infty}.$

Analysis of (S, R)

We denote by $R_0 = \beta/\gamma$ which is the famous Basic Reproduction Number of the disease. As $S(t) + I(t) + R(t) = N$ , for all $t>0$, we can eliminate $I$ and obtain a system of $S, R$ only.

From the exercise, we can let $t\to \infty$ and the deriavative become zero as the functions apporache to constants. Therefore the right hand side of $\eqref{SRN}$ is also zero. We then obtain the relation $S_{\infty}+ R_{\infty} = N$ and thus $I_{\infty} = \lim_{t\to \infty} I(t) = 0$. We thus proved obersvation 1 and 2.

Let's try to figure out an equation for $S, R$. The first equation in (12) divided by the second one imples an ODE

with initial condition $S(0), R(0) = 0$. It can be easily solved to get

Substitute $S$ in terms of $R$ in the second equation of $\eqref{SRN}$, we get a nonlinear ODE of $R$

Unfortunately I can't solve this onlinear ODE.

But we can get an algebraic equation of the limit. Let $t\to \infty$ and use $R_{\infty} = N - S_{\infty}$, we obtain from (13) the relation

whose solution can be easily obtained, e.g., by Newton's method. In application, $I_0 \ll N$ ( e.g. 1 patient over million people) and $S(0)\approx N$. With an abuse of notation, we set $S_{\infty} \leftarrow S_{\infty}/N$, we can solve the equation

to get an accurate estimate of the population ratio which remains non-infected in the end.

Analysis of I

We turn to the estimate of infected number $I(t)$. We rewrite the equation as

We can solve $I$ in terms of $S$ as

Again we only obtain the relation between $I(t)$ and $S(t)$ not analytic formulation of the solution.

Remember $S(0)\approx N$ and thus $S/N \approx 1$ in the beginning and $S/N < 1$ always true.

The behavior of $I$ depends on crucially $R_0$

1. If $R_0 \leq 1$, then $I' < 0$ which means the number of infected dies off exponetially. Recall that $R_0 = \beta/\gamma$. It is less than one means the recover rate is greater than the infection rate. The disease will disappear. Indeed in this case, the disease is not considered as epidemic.

2. We then consider the epidemic case $R_0 >1$. From the first equation of SIR model, $S$ is decreasing.

   • when $\frac{S}{N} > \frac{1}{R_0}$, then $I' > 0$, i.e., $I$ is increasing.
   • the maximum infected number is acheived at the time when $S(t)/N = 1/R_0$.
   • after that, $I' < 0$, i.e., $I$ is decreasing to zero.

We explore more on these 3 stages. We first plot $I$ in a different figure.

1. In the early stage, $S/N \approx 1$, by $\eqref{IS}$, the infected number is

   $$I(t) \approx I_0\exp(\gamma (R_0 - 1)t) = I_0 \exp ((\beta - \gamma)t).$$

   Namely we obersve the exponential increase of infectious in the early stage. The rate is mainly determined by $\beta$ offside by $\gamma$.

2. When $S \leq N/R_0$​ , the population starts having enough of what is known as Herd immunity for the disease to be unable to spread further.

3. We calculate the peak of $I$ as follows. Use the relation (13) for $R$ and $S$ , and use $S(0)\approx N$, we get, when $\frac{S(t)}{N} = \frac{1}{R_0}, \frac{R(t)}{N} = \frac{\ln R_0}{R_0}$ , combine with the conservation of population $S(t) + I(t) + R(t) = N$ , we get

   $$\frac{I_{\max}}{N} = 1 - \frac{1+\ln R_0}{R_0}.$$

   With this formulation, it is easy to see $I_{\max}$ is an increasing function of $R_0$. Flatten the curve means reducing the peak $I_{\max}$ below the health care system capacity.

   Exercise Estimate the peak time $t^*$.

   

    

Effect of Social distancing

Ways to reduce $R_0 = \beta/\gamma$

• increae the recovery rate $\gamma$ by medicine
• reduce the infection rate $\beta$ by vacinne or by social distancing.

Parameters for the figure: $N = 1000, I_0 = 10, \beta = 0.1, \gamma = 0.05$. When $t>50$, we set $\beta = 0.01$ and show the results by the dashed curves. The crticial parameter $R_0$ is dropped from $2$ to $0.2$. We see immediately the infectious starts to decay exponentially and then the majority remains non-infected almost immediately (exponentially fast).

For coronavirus, the change is not that dramatically. This is due to the time delay for those who exposed to the infected but get infected after 7 - 14 days.

So we need to add one more compartment: Exposed and discuss SEIR model in the next section.

SEIR model

We introduce one more compartment $E$ and one more parameter $\alpha$ and list them below.

• $S(t)$ ratio of the susceptible population;
• $E(t)$ ratio of the exposed population;
• $I(t)$ ratio of the infected population;
• $R(t)$ ratio of the recovered population

Parameters $(\alpha, \beta, \gamma)$.

• $\alpha$ is the inverse of the incubation period
• $\beta$ is the contact rate
• $\gamma$ is the recovery rate.

The basic reproduction number is still defined by $R_0 = \beta/\gamma$ .

Model

We add $E$ as a transition between $S$ and $I$ and the system reads as

We use $\dot{S}$ to denote the derivative with respect to $t$. More importantly, we normlize all functions by dividing $N$ and thus it represents the ratio to the whole pouplation. Conservation of population becomes $S(t) + E(t) + I(t) + R(t) = 1$ for $t>0$ which can be preserved by chosing initial condition satisfying $S(0) + E(0) + I(0) + R(0) = 1$. It is reasonable to chose $I(0) =R(0) = 0$ and $E(0) \ll 1$ and thus $S(0)\approx 1$, e.g., $S(0) = 0.9999$.

Numerical solution

We write out the forward Euler's method with $\Delta t = 1$.

which gives a better explanation of the model. For example, the last one means tomorrow's recovered people is equal to today's add the recovered ratio $\gamma$ multiply the number of infected people $I$.

We use $\gamma = 0.8$ as according many report $80\%$ of infected people will recovery without hosptilization. The rest $20\%$ will go to hosptial. We use $R_0 = 3.5$ which is reasonable according several studies. So $\beta = 2.8$. We use the averaged incubation period as $5$ days which implies $\alpha = 0.2$.

Although the peak of $I$ is relatively small (thanks to the large recovery rate), eventually over $95\%$ will be infected due to the large $R_0$.

Again we consider strick rules to reduce $R_0$ to $R_0/10$. The first figure is control time starting from $30$ days and the second one is $25$ days. The change is dramatical.

Only $5$ days earlier, but the recovered popluation is reduced from almost $60\%$ to $25\%$.

Perturbated Analysis

We discuss perturbated analysis for solving nonlinear ODEs. We use normalized variables $S, I, R \in [0,1]$ and the nonlinear ODE we obtained is

with initial condition $R(0) = 0$ and $0 < 1 - S(0) \ll 1$.

Let $f(x) = 1 - S(0)e^{-R_0 x} - x$. Then solution to $\eqref{nonlinearR}$ can be written as

The problem is we do not the antiderative of $1/f(x)$ for such a nonlinear $f(x)$.

We shall use polynomial approximations to $f(x)$ at different point $x$ to find approximation of $\eqref{solution}$ . First of all, we use approximation $S(0) = 1$ to simplify $f(x) = 1 - e^{R_0 x} - x$. Then we build the condition $\bar R = R(\bar t)$ into the solution in the form, i.e., integrate from $(\bar t, \bar R)$,

Linear approximation

Consider $f(x) = a x + b$. We get the solution

Quadratic approximation

Consider $f(x) = ax^2 + b x + c = a (x - \lambda_1)(x - \lambda_2)$ where $\lambda_1\leq \lambda_2$ are two real roots of $f(x)$. Use the factorization

we get

and

which is a logistic function transition from $\lambda_1$ to $\lambda_2$ and the critical point, i.e., $R't^*) = 0$ is given by

We apply the formula to several points and get approximation at different stages.

Initial stage

We apply the formula to $\bar t=0, \bar R = 0$. For the linear approximation, we obtain

Use $e^{-x}\approx 1 - x + \frac{1}{2}x^2$ to get a quadratic ODE at $x=0$. The corresponding setting is

Ending stage

We apply the formula at $\bar t = t_{\infty}, \bar R = R_{\infty}$ . For the linear approximation,

For the quadratic approximation, we obtain

Turning point stage

We apply the formula at $\bar t = t^*, \bar R = R^*$ , where $t^*$ is the peak time where $I$ achieves the maximum, i.e. $I(t^*) = I_{\max}$. From the analysis of SI model, we know

Reference

1. Epidemic Modeling 101: Or why your CoVID-19 exponential fits are wrong https://medium.com/data-for-science/epidemic-modeling-101-or-why-your-covid19-exponential-fits-are-wrong-97aa50c55f8
2. 知乎：关于传染病的数学模型 https://www.zhihu.com/question/367466399