Approximated solutions of SEIR models

Question

Numerical solutions of the SEIR equations (describing the spreading of an epidemic disease) – or variations thereof –

• $\dot{S} = - N$

• $\dot{E} = + N - E/\lambda$

• $\dot{I} = + E/\lambda - I/\delta$

• $\dot{R} = + I/\delta$

with

• $N = \beta I S / M$ = number of newly infected individuals

• $\beta = $ infection rate

• $\lambda = $ latency period

• $\delta = $ duration of infectiosity

• $M = S + E + I + R = $ size of the population

yield characteristic and almost symmetric peaks for the function $I(t)$ of numbers of infectious individuals. So $I(t)$ can – by a rough guess – be approximated by a Gauss curve

$$\widetilde{I}(t) = I_0\ \operatorname{exp}\Big({-\big((t-t_0)/\sigma\big)^2}\Big)$$

with $I_0$ the maximal value of $I(t)$, $I(t_0) = I_0$, and $\sigma$ such that $\widetilde{I}(0) = I(0) = 1$, i.e.

$$\sigma = t_0\ /\ \sqrt{\text{ln} I_0}$$

For different values of $\delta$, the reproduction number $R_0 = \beta\cdot\delta$, and a fixed value $\lambda = 2$ we find:

It turns out that an exponent $\sqrt{2}$ instead of $2$ yields better results, i.e.

$$\widetilde{I}(t) = I_0\ \operatorname{exp}\Big({-\big(|t-t_0|/\sigma\big)^{\sqrt{2}}}\Big)$$

My question is fourfold:

1. Why is a Gauss-like curve a good approximation at all? That means: Why is $I(t)$ so symmetric?

2. By which considerations could one come up with the exponent $\approx \sqrt{2}$?

3. By which considerations can the asymmetry of the numerical solution $I(t)$ be understood which becomes apparent when comparing it with the symmetric approximation $\tilde{I}(t)$?

4. Has anyone an idea how $I_0$ and $t_0$ look like as functions of $\beta,\lambda,\delta,M$?

---

Just to give another view on the tables above, find here all curves overlayed:

Answer 1

Is your function F(x) known for giving good fits in many contexts? How can these be characterized?

This is too long for a comment but I'd like you to check if the fit is to your satisfaction before I elaborate. I prefer to write everything in the numerator, so my equations will be $$ \dot S=-\beta IS, \dot E=\beta IS-\lambda E, \dot I=\lambda E-\delta I\,. $$ Suppose that $I_0$ is the maximum of $I$ attained at the moment $0$ (just shift otherwise).
Then the equations I'm using (I hope I'm copyng them right) are $$ 2a^2(\beta I_0+\lambda+\delta-\mu)=\lambda\delta \beta I_0 \\ 6a^2=(\lambda+\delta-\mu)(\beta I_0-\mu) $$ Once you have solved those for $a,\mu>0$ (assume that $I_0$ is known for the moment and you just want a fitting curve rather than an independent derivation for everything), let $a_\pm=\sqrt{a^2+\frac{\mu^2}4}\mp \frac\mu 2$ (so $a_->a_+$), define $$ F_{a,\mu}(t)=I_0\left(\frac{a_-\exp(a_+t)+a_+\exp(-a_-t)}{a_-+a_+}\right)^{-2} $$ and compare it to $I(t)$. If you like the fit, we can discuss where all that nonsense came from and how to write the full system where $I_0$ will be solved for, not given. If not, I'll stop here, so let me know what you think.

The equations are algebraic of third degree, so, unless you are a big fan of Cardano's formulae, you'll have to solve them numerically. That's not hard (almost any decent iteration scheme works). The approximation is pretty good in most cases, IMHO, but it has its limitations so one can find regimes where it breaks though those are usually rather extreme. Enjoy! :-)

Two pictures, as promised. The black curve is the true trajectory, the red one is the computed trajectory (note that the height of the peak is also computed: I finally found a good third equation, so I played it honestly and didn't try to tweak the parameters beyond what my linearized equations gave directly), the green line is the best symmetric approximation you can hope for (half sum of the true trajectory and its reflection around the peak). I believe that the red line is better even without any tweaking and that the precision with which the maximum is determined is also fairly decent, but you can judge by yourself :-).

Answer 2

From an article that user @Gro-Tsen refers to I learned - and give here as a partial answer - that for the case of a vanishing latency period $\lambda = 0$, i.e. for the classical SIR model

• $\dot{S} = -N$
• $\dot{I} = +N - I/\delta$
• $\dot{R} = +I/\delta$

there is a closed formula for $I_{max}$ (i.e. the maximal value of $I(t)$) as a function of $\beta$, $\delta$, and $M$ namely

$$ I_{max} = \frac{R_0 - \log R_0 - 1}{R_0} \cdot M$$

with $R_0 = \beta\cdot\delta$. This is quite nice.