Since we are currently fighting the Corona epidemy I thought it might be a good idea to simulate the spreading of such a disease in Euler Math Toolbox. The numbers here are fictitious and can only deliver qualitative results.

The mathematical model I took is quite simple.

• The disease has n1=7 days of incubation. During this time the infected person can infect others.
• The disease has n2=7 days of severe sickness. I have assumed that there are no further infections during that time.
• We start with the a-immune population plus 10 infected cases.
• The chance for infection depends on the probability to contact an infected person during the incubation time and the number of contacts each day that could lead to an infection. I assume that a person has nc=0.2 to nc=0.3 contacts each day which could lead to an infection.
• In the end, the population is either in the non-immune or the immune group.

So the total population is mapped into a vector x of groups as

• x[1] : non-immune
• x[2] to x[1+n1] : incubation
• x[2+n1] to x[1+n1+n2] : sick
• x[2+n1+n2] : immune

The formula for the number of new infections in x[2] (first day of incubation) is

\(x[2] = x[1] \left(1-(1-\frac{S}{P})^{n_c}\right)\)

where S is the number of incubated people and P is the total population. So 1-S/P is the chance that the contacted person is not infected, and thus the factor computes the chance for infection. Note that this factor is very small at the start since S/P is very small. Yet the exponential growth does its job over time.

Here is a typical plot of the epidemy with nc=0.3. This is equivalent to contacting 9 other people in a month in a way that can spread the disease. We see that almost 90% of the population has to go through the disease in this case.

Here is another plot with nc=0.25, equivalent to contacting 7.5 other people in a month.

That reduces the number of cases a lot and limits the peak number of sick persons from 37% to 22%. For the medical system, this is an enormous relief.

The code of Euler Math Toolbox to do these plots is listed below. You can cut and paste it into EMT to try yourself.

>n1=7; // 7 days of incubation time
>n2=7; // 7 days of sickness
>pm=0.01; // mortility rate
>pop=50000000; // population
>nc=0.3; // number of infectious contacts
>function oneday (x,n1,n2,pm,pop,nc) ...
$xn=x;
$xn[2]=floor(xn[1]*(1-(1-sum(x[2:1+n1])/pop)^nc)); // new infections
$xn[1]=x[1]-xn[2]; // remaining uninfected
$xn[3:2+n1+n2]=x[2:1+n1+n2]; // progress one day 
$xn[-1]=xn[-1]+x[-1];
$return xn;
$endfunction
>x=zeros(2+n1+n2); x[1]=pop; x[2]=10; // start with 10 cases
>function sim (x,n1,n2,pm,pop,nc) ...
$X=x;
$repeat
$   xn=oneday(x,n1,n2,pm,pop,nc);
$   X=X_xn;   
$   until sum(xn[2:1+n1+n2])==0;
$   x=xn;
$end;
$return X;
$endfunction
>X=sim(x,n1,n2,pm,pop,nc);
>Xs=X[,1]|sum(X[,2:2+n1])|sum(X[,2+n1:1+n1+n2])|X[,-1];
>plot2d(Xs'/pop,color=[green,orange,red,blue]); ...
>labelbox(["immune-","infect+","sick","immune+"], ...
>  colors=[green,orange,red,blue],x=0.4,y=0.2,w=0.35):

I do have the same code in Python.

import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt

n1=7 # 7 days of incubation
n2=7 # 7 days of sickness
pm=0.01 # mortility rate
pop=50000000 # population
nc=0.3 # number of infectious contacts

def oneday (x,n1,n2,pm,pop,nc):
    xn=x.copy()
    xn[1]=np.floor(x[0]*(1-(1-sum(x[1:1+n1])/pop)**nc)) # new infections
    xn[0]=x[0]-xn[1] # remaining uninfected
    xn[2:2+n1+n2]=x[1:1+n1+n2] # progress one day 
    xn[-1]=xn[-1]+x[-1]
    return xn

x=np.zeros(2+n1+n2)
x[0]=pop
x[1]=10 # start with 10 cases

def sim (x,n1,n2,pm,pop,nc):
    X=x
    while True:
        xn=oneday(x,n1,n2,pm,pop,nc)
        X=np.vstack((X,xn))
        if np.sum(xn[2:1+n1+n2])==0:
            break
        x=xn;
    return X;

X=sim(x,n1,n2,pm,pop,nc)

x1 = X[:,0]
x2 = np.sum(X[:,1:1+n1],axis=1)
x3 = np.sum(X[:,1+n1:n1+n2],axis=1)
x4 = X[:,-1]

fig,ax = plt.subplots(figsize=(14,8))
ax.grid(True)
ax.plot(x1,label='immune-')
ax.plot(x2,label='infect+')
ax.plot(x3,label='sick')
ax.plot(x4,label='immune+')

plt.legend();

We can also understand the problem a bit better if we study the number of new infections in x[2], the first day of incubation. It starts with our chosen number of initial infections. But those 10 induce 3 new infections and the 13 new take another 3. Then we have 16 and they infect 4, and so on exponentially.

If the chance S/P of an infected contact is small we have

\(\left(1-(1-\frac{S}{P})^{n_c}\right) \sim \frac{n_c S}{P}\)

Thus

\(x[2] \sim n_c S\)

exactly as observed.