Perturbation Analysis of SIR Model

SIR Model

Model

We consider SIR model which divides the population into three compartments:

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

All three variables are normlized by dividing the total population $N$ so that their values are in $[0,1]$. Introduce two parameters:

• $\beta$ : the infection rate;
• $\gamma$ : the recovery rate.

The model is

with initial condition $S(0) = S_0 \approx 1, I(0) = I_0$, and $I_0 + S_0 + R(0) = 1$. Note that the notation $R_0 = \beta/\gamma > 1$ is the famous Basic Reproduction Number of the disease not the initial value of $R$.

Explanation: the product $IS$ is the number of susceptible person will encounter an infected one. Multiply the rate of infection $\beta$ is the new infected population. $\gamma I$ is the rate of recovered population.

Simply add these three equations, we have $\frac{{\rm d}}{{\rm d} t}(S(t) + I(t) + R(t)) = 0$ and thus $S(t) + I(t) + R(t) = 1$ for all $t$. This is known as the conservation of population. No death or birth is included in the model.

Analysis of (S, R)

Since $S(t) + I(t) + R(t) = 1$ , for all $t>0$, we can eliminate $I$ and obtain a system of $S, R$ only.

The first equation of $\eqref{SIR}$ implies $S$ is decreasing and the last one implies $R$ is increasing. As they are bounded in $[0,1]$, their limit exist as $t\to \infty$ and denote by $S_{\infty}, R_{\infty}$.

We then let $t\to \infty$ and the deriavatives on the left hand side of $\eqref{SR}$ become zero as the functions apporache to constants. Therefore the right hand side of $\eqref{SR}$ should also be zero. We then obtain the relation $S_{\infty}+ R_{\infty} = 1$ and consequently $I_{\infty} = \lim_{t\to \infty} I(t) = 0$.

Let's try to solve $\eqref{SR}$. Divide the first equation in $\eqref{SR}$ by the second one imples an ODE

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

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

Let $t\to \infty$ and use $R_{\infty} = 1 - S_{\infty}$, we obtain from $\eqref{SRrelation}$ the relation

whose solution can be easily obtained by, e.g., Newton's method. In real application, $I_0 \ll 1$ ( e.g. 1 patient over million people) and we can replace $S_0$ by $1$ and solve the equation

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

Analysis of I

We rewrite the equation for $I$ as

As $S$ is decreasing and in $(0,1)$ and for the epidemic model $R_0>1$ , we can divide into three stages according to the sign of $S - 1/R_0$ .

1. Early stage: when $S > \frac{1}{R_0}$, then $I' > 0$, i.e., $I$ is increasing.
2. Peak stage: the maximum infected number is acheived at the time $t^*$ when $S(t^*) = 1/R_0$.
3. Ending stage: after the peak time $I' < 0$, then $I$ is decreasing to zero.

Perturbation Analysis

We apply perturbation analysis for solving nonlinear ODEs. We use normalized variables $S, I, R \in [0,1]$ and replace $S_0$ by 1. The nonlinear ODE we try to solve is

with an apporiate initial condition . Let $f(x) = 1 - e^{-R_0 x} - x$. From the third equation of SIR model $R'=\gamma I$, we have the relation $I(t) =f(R(t))$.

Then solution to $\eqref{nonlinearR}$ can be written as

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

What can we solve? We can solve when $f$ is simple polynomials such as linear or quadratic. So we can use polynomial approximations to $f(x)$ at different point $x$ to find approximations of $\eqref{solution}$ . More specifically, given time $\bar t$ and a value $\bar R\in (0,1)$, we build the condition $\bar R = R(\bar t)$ into the solution by the following form

Note that we have to assume the interval $[\bar R, R]$ does not contains the root of $f(x)$. Otherwise the integeral on the left could be divergence.

Linear approximation

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

Note that when $b=0$, we cannot integrate $1/x$ from $0$. That is $\bar R > 0$ and cannot be zero. Using $R'=\gamma I$, we get

So we can use $\bar I$ to determine $\bar R$.

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)$. Here we consider only the case the quadratic function has two real roots. To be well defined, we assume the interval $(\bar R, R) \subset [\lambda_1, \lambda_2]$ so that the roots of $f(x)$ are not in $(\bar R, R)$.

Use the factorization

and integrate $\eqref{solution2}$, we get

where rate $\quad \rho = a\gamma (\lambda_1 - \lambda_2)$. We can solve $\eqref{Rfirst}$ to get

with parameters

which is a logistic function transition from $\lambda_1$ to $\lambda_2$ as the limits

Its critical point, i.e., $R_Q'(t^*) = 0$ is given by

Using the relation $R' = \gamma I$, we get the formulae on $I_Q$

Taylor expansion

We list the function and its derivatives

The Taylor series at a certain point $\bar x$ is

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

Initial stage

We apply the formula to $\bar t=0, \bar x = 0$. Then

For the linear approximation $a=R_0-1$ and $b=0$, we obtain

Namely in the initial stage, it is like the solution of ODE $y' = \beta y - \gamma y$. The solution grows with rate $\beta$ and decays with rate $\gamma$.

As mentioned early, $\bar R >0$ and cannot be zero. We can determine $\bar R$ by looking at $I(t) = R'(t)/\gamma$ for which the initial condition $I_0$ is given. From

we get $\bar R = I_0/(R_0-1)$ which is small but positive. That is $\bar R$ is a small perturbation of $0$ but not zero.

The quadratic approximation of $f$ near $x=0$ is $(R_0 - 1) x - \frac{1}{2}R_0^2x^2 = x(R_0 - 1 - \frac{1}{2}R_0^2x)$ which implies

Instead of pluging into the formula to simplify, we leave the general formula as it can be computed using values of $a, \lambda_1,\lambda_2$. Note that the rate $\rho = \beta -\gamma$ is the same as the one in the linear approximation.

We compare the analytic approximation with the numerical solution. In these figures, blue curves are numerical solution obtained by ode45 in MATLAB and red is the linear approximation and black is the quadratic one.

In the early stage, e.g, $0\leq t \leq 30$, both linear and quadratic approximations are very accurate. But after that, the linear approximation keeps exponentially increasing with the same rate while the quadratic one will slow down and approaches to a limit (for $R$) or zero (for $I$) due to the behavior of logistic functions.

We can use formula $\eqref{tbartstar}$ to predict the peak time but be aware that could be optimistic. The infected number could still increase for a while.

Ending stage

We compute derivatives of $f(x)$ at $\bar x = R_{\infty}$

The function value $f(R_{\infty}) = I(\infty) = 0$ . On the other hand, by the formulation of $f(x)$, we obtain the relation $\exp(-R_0 R_{\infty}) = 1-R_{\infty} = S_{\infty}$ and thus the derivatives can be simplified to

For the linear approximation, we have $a = R_0S_{\infty} - 1, \, b/a = - R_{\infty}$. The linear approximation of $R$ is

and linear approximation of $I$ is

with rate $\rho = \gamma - \beta S_{\infty}$. Recall that the increasing rate in the early stage is $\beta - \gamma$ which is usally bigger than the decay rate in the end.

If we chose $\bar R = R_{\infty}$, we get constant approximations which are not interesting. Instead we chose $\bar t, \bar R = R(\bar t)$ when we observe the function $I$ starts to decay for a while but close to the flat tail. The linear approximation tells us although the decay is not as fast as the pandemic starts, it will decay to zero exponentially.

For the quadratic approximation, we have the Taylor expansion $f'(R_{\infty})(x - R_{\infty}) + \frac{1}{2} f''(R_{\infty})(x-R_{\infty})^2 = \frac{1}{2} f''(R_{\infty})(x-R_{\infty})\left [2f'(R_{\infty})/f''(R_{\infty})+ x - R_{\infty}\right ].$

Therefore

Again we plug into the general formulation and compare with the numerical solution.

Again the quadratic one provides a much better approximation In the ending stage and not trustable in the initial stage.

Turning point stage

We shall apply the formula at $\bar t = t^*, \bar x = \bar R = R^*$ , where $t^*$ is the peak time when $I$ achieves the maximum. Denote by $S^* = S(t^*), I^* = I(t^*), R^* = R(t^*).$

From the equation for $I$,

by solving $I'(t^*) = 0$, we know

Use relation $I = f(R)$, we can easily get

Therefore the linear approximation of $I$ is the constant approximation $I_L = I^*$ and

The quadratic approximation is given by

from which we get

By the symmetry, we have simplified formulae

and

The comparison to the numerical solution can be found in the figures below. The approximation $I_Q$ is surprisingly good for all time $t$. It is worthing to point out the shape of the function is determined by the single parameter $R_0$. To determine the location, the peak time $t^*$ is also needed.

Global approximation

We can glue quadratic approximation in different stages to get a piecewise logistic function approximation of the true solution of SIR model although we cannot solve it analytically.

We can assume the peak time $t^*$ is observable from the data. We chose a small interval near $t^*$ like $[t^* - t_1, t^* + t_2]$. Denote the quadratic approximation in the early, peak, and ending stage as $R_{Q}^1,R_{Q}^2,R_{Q}^3$, respectively. We will use the piecewise approximation

Similar definition is applied to obtain $I_g$.

We plot the approximation $R_g, I_g$ by black dash lines which are almost not distinguishable from the numerical solution.