Simple model of infection

Each individual infects $c$ new people each day.

$I_{n+1} = I_n + (c I_n) = I_n(1 + c)$

So

$I_{n+1} = \lambda I_n$, where $\lambda := 1 + c$ -- growth rate

and $I_0 = 1$

In [3]:
I0 = 1
Out[3]:
1
In [4]:
I_0 = 1
Out[4]:
1
In [6]:
I₀ = 1   # \_0 <TAB>
Out[6]:
1
In [7]:
c = 0.01  # average no. of people that each individual infects on each day
Out[7]:
0.01
In [8]:
λ = 1 + c  # \lambda <TAB>
Out[8]:
1.01
In [10]:
λ = 1 + c;

; suppresses output

In [11]:
I_1 = λ * I_0
Out[11]:
1.01
In [12]:
I_2 = λ * I_1
Out[12]:
1.0201
In [13]:
I_3 = λ * I_2
Out[13]:
1.030301

Arrays

In [14]:
T = 10  # final time 

I = zeros(T)
Out[14]:
10-element Array{Float64,1}:
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
In [15]:
I
Out[15]:
10-element Array{Float64,1}:
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
In [16]:
I[1]
Out[16]:
0.0
In [17]:
I[1:2]
Out[17]:
2-element Array{Float64,1}:
 0.0
 0.0
In [19]:
I = zeros(Int64, T)
Out[19]:
10-element Array{Int64,1}:
 0
 0
 0
 0
 0
 0
 0
 0
 0
 0
In [20]:
I = zeros(T)
Out[20]:
10-element Array{Float64,1}:
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
In [21]:
I[1] = I_0   # Could use OffsetArrays.jl -- enables arbitrary indexing
Out[21]:
1
In [25]:
for n in 1:T-1
    I[n+1] = λ * I[n]
    @show n, I[n]
end

# for loops do not return anything so running a for loop does not output
(n, I[n]) = (1, 1.0)
(n, I[n]) = (2, 1.01)
(n, I[n]) = (3, 1.0201)
(n, I[n]) = (4, 1.030301)
(n, I[n]) = (5, 1.04060401)
(n, I[n]) = (6, 1.0510100501)
(n, I[n]) = (7, 1.061520150601)
(n, I[n]) = (8, 1.0721353521070098)
(n, I[n]) = (9, 1.08285670562808)
In [26]:
I
Out[26]:
10-element Array{Float64,1}:
 1.0
 1.01
 1.0201
 1.030301
 1.04060401
 1.0510100501
 1.061520150601
 1.0721353521070098
 1.08285670562808
 1.0936852726843609
In [27]:
using Plots
┌ Info: Precompiling Plots [91a5bcdd-55d7-5caf-9e0b-520d859cae80]
└ @ Base loading.jl:1260
In [29]:
plot(I, m=:o, label="I[n]", legend=:topleft)
Out[29]:
2 4 6 8 10 1.00 1.02 1.04 1.06 1.08 I[n]
In [30]:
T = 20

I = zeros(I)
I[1] = I_0

for n in 1:T-1
    I[n+1] = λ * I[n]
    @show n, I[n]
end

plot(I)
(n, I[n]) = (1, 1.0)
(n, I[n]) = (2, 1.01)
(n, I[n]) = (3, 1.0201)
(n, I[n]) = (4, 1.030301)
(n, I[n]) = (5, 1.04060401)
(n, I[n]) = (6, 1.0510100501)
(n, I[n]) = (7, 1.061520150601)
(n, I[n]) = (8, 1.0721353521070098)
(n, I[n]) = (9, 1.08285670562808)
BoundsError: attempt to access 10-element Array{Float64,1} at index [11]

Stacktrace:
 [1] setindex!(::Array{Float64,1}, ::Float64, ::Int64) at ./array.jl:825
 [2] top-level scope at ./In[30]:6

Functions

In [31]:
function run_infection(T=20)  # default

    I = zeros(T)
    I[1] = I_0

    for n in 1:T-1
        I[n+1] = λ * I[n]
    end

    return I
end
Out[31]:
run_infection (generic function with 2 methods)
In [32]:
methods(run_infection)
Out[32]:
# 2 methods for generic function run_infection:
  • run_infection() in Main at In[31]:3
  • run_infection(T) in Main at In[31]:3
In [33]:
run_infection(10)
Out[33]:
10-element Array{Float64,1}:
 1.0
 1.01
 1.0201
 1.030301
 1.04060401
 1.0510100501
 1.061520150601
 1.0721353521070098
 1.08285670562808
 1.0936852726843609
In [34]:
I_result = run_infection(10);
In [35]:
I_result
Out[35]:
10-element Array{Float64,1}:
 1.0
 1.01
 1.0201
 1.030301
 1.04060401
 1.0510100501
 1.061520150601
 1.0721353521070098
 1.08285670562808
 1.0936852726843609

Always separate data generation from plotting!

In [38]:
I_result = run_infection(20)

plot(I_result, m=:o)
Out[38]:
5 10 15 20 1.00 1.05 1.10 1.15 1.20 y1
In [39]:
using Interact

Unable to load WebIO. Please make sure WebIO works for your Jupyter client. For troubleshooting, please see the WebIO/IJulia documentation.

In [44]:
end_T = 1000
@manipulate for T in slider(1:end_T, value=1)
    I_result = run_infection(T)

    plot(I_result, m=:o)
    
    xlims!(0, end_T)
    ylims!(0, 10)
end
Out[44]:

Exponential growth

In [46]:
plot(run_infection(1000))
Out[46]:
0 250 500 750 1000 0 5.0×10 3 1.0×10 4 1.5×10 4 2.0×10 4 y1
In [47]:
plot(run_infection(1000), yscale=:log10)
Out[47]:
0 250 500 750 1000 10 0 10 1 10 2 10 3 10 4 y1
In [50]:
I_result = run_infection(1000)

plot(log10.(I_result))   # log of each element in I_result
ylabel!("log(I_n)")
Out[50]:
0 250 500 750 1000 0 1 2 3 4 log(I_n) y1

$\log(I_n) = a n + b$

Take exponentials of both sides:

$I_n = \exp(an + b) = C e^{a n}$

Solve $I_{n+1} = \lambda I_n$:

$I_n = \lambda^n I_0$

$\log(\lambda^n) = n \log(\lambda)$

Logistic growth

Exp growth is unrealistic: Assumes that there are always more people to infect -- wrong (finite population).

$I_{n+1} = I_n + c I_n$

Each individual will be in contact with a fraction $\alpha$ of the population. At each contact there will a probability $p$ that you infect each person.

Number of contacts = $\alpha N$

So $c = p \alpha N$

Original model:

$I_{n+1} = I_n + (p \alpha N) I_n$ -- good approximation when almost everybody still susceptible

New model:

Can only infect uninfected people!

$I_{n+1} = I_n + (p \alpha S_n) I_n$

$I_{n+1} = I_n + [p \alpha (N - I_n)] I_n = f(I_n)$

$I_{n+1} = I_n + \beta(I_n, S_n) I_n = f(I_n)$

In [56]:
p = 0.02
α = 0.01
N = 1000

β(I, S) = p * α * S

function run_infection(T=20)  # default

    I = zeros(T)
    I[1] = I_0

    for n in 1:T-1
        I[n+1] = I[n] + β(I[n], N - I[n]) * I[n]
    end

    return I
end
Out[56]:
run_infection (generic function with 2 methods)
In [57]:
I = run_infection(20)
Out[57]:
20-element Array{Float64,1}:
  1.0
  1.1998
  1.439472095992
  1.7269520992073721
  2.071746046338255
  2.4852368292698026
  2.9810489147042514
  3.57548136711873
  4.288020827141153
  5.141947568046585
  6.165049156697394
  7.390457421815974
  8.857625133998434
 10.613458656195233
 12.713621286504926
 15.224018310562577
 18.222467825971023
 21.80054972443152
 26.065606875660343
 31.14284507843312
In [58]:
plot(I, m=:o)
Out[58]:
5 10 15 20 5 10 15 20 25 30 y1
In [59]:
plot(I, m=:o, yscale=:log10)
Out[59]:
5 10 15 20 10 0.0 10 0.5 10 1.0 10 1.5 y1
In [61]:
I = run_infection(100)
plot(I, m=:o, yscale=:log10)
Out[61]:
0 25 50 75 100 10 0 10 1 10 2 10 3 y1
In [64]:
p = 0.01
α = 0.1
N = 100

I = run_infection(100)
plot(I, m=:o)
Out[64]:
0 25 50 75 100 0 25 50 75 100 y1

Sigmoid shape -- S-shaped

In [65]:
plot(I, m=:o, yscale=:log10)
Out[65]:
0 25 50 75 100 10 0.0 10 0.5 10 1.0 10 1.5 10 2.0 y1

"Logistic-type growth"

$I_{n+1} = \lambda I_n$

More realism: Heterogeneity of individuals or groups

Idea: Instead of modelling population globally, model each individual

"Patch model": Local patches where population is well mixed

Network model: Links between nodes

Exponential growth:

$I_{n+1} = \lambda_n I_n$ -- growth rate changes in time

In [66]:
rand()
Out[66]:
0.2767437763441829
In [67]:
rand()
Out[67]:
0.7089008345855983
In [68]:
rand()
Out[68]:
0.3789568897206148

www.random.org

In [69]:
randn()
Out[69]:
-0.33523291258515997
In [ ]:

In [71]:
p = 0.02
α = 0.01
N = 1000



function run_infection(c_average=1.1, T=20)  # default

    I = zeros(T)
    I[1] = I_0

    for n in 1:T-1
        c = c_average + 0.1*randn()
        I[n+1] = I[n] + c * I[n]
    end

    return I
end
Out[71]:
run_infection (generic function with 3 methods)
In [73]:
c_average = 1.1

cs = [c_average + 0.1*randn() for i in 1:100]
Out[73]:
100-element Array{Float64,1}:
 1.1794683814069813
 1.0510617699472489
 1.0799811965453343
 1.3160413008055867
 1.0655215426073297
 1.014913328568307
 1.2142141455660913
 0.9952511944053644
 1.1118228511102424
 1.2274911767386665
 1.097194489799497
 1.156157408743076
 1.0138054185809025
 ⋮
 1.252209922391299
 1.0305288124096037
 1.1816089517166688
 0.9848866930013738
 1.175520029029598
 1.1695318311722314
 1.1237453859981914
 1.1925550201637143
 0.9331362892813453
 1.0848564135882912
 1.184898469695791
 1.0756350756110384
In [74]:
scatter(cs)
Out[74]:
0 25 50 75 100 0.8 0.9 1.0 1.1 1.2 1.3 y1
In [ ]: