# 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]

In [29]:
plot(I, m=:o, label="I[n]", legend=:topleft)

Out[29]:
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]:
In [39]:
using Interact


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]:
In [47]:
plot(run_infection(1000), yscale=:log10)

Out[47]:
In [50]:
I_result = run_infection(1000)

plot(log10.(I_result))   # log of each element in I_result
ylabel!("log(I_n)")

Out[50]:

$\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]:
In [59]:
plot(I, m=:o, yscale=:log10)

Out[59]:
In [61]:
I = run_infection(100)
plot(I, m=:o, yscale=:log10)

Out[61]:
In [64]:
p = 0.01
α = 0.1
N = 100

I = run_infection(100)
plot(I, m=:o)

Out[64]:

Sigmoid shape -- S-shaped

In [65]:
plot(I, m=:o, yscale=:log10)

Out[65]:

"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

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]:
In [ ]: