# 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 :
I0 = 1

Out:
1
In :
I_0 = 1

Out:
1
In :
I₀ = 1   # \_0 <TAB>

Out:
1
In :
c = 0.01  # average no. of people that each individual infects on each day

Out:
0.01
In :
λ = 1 + c  # \lambda <TAB>

Out:
1.01
In :
λ = 1 + c;


; suppresses output

In :
I_1 = λ * I_0

Out:
1.01
In :
I_2 = λ * I_1

Out:
1.0201
In :
I_3 = λ * I_2

Out:
1.030301

## Arrays¶

In :
T = 10  # final time

I = zeros(T)

Out:
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 :
I

Out:
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 :
I

Out:
0.0
In :
I[1:2]

Out:
2-element Array{Float64,1}:
0.0
0.0
In :
I = zeros(Int64, T)

Out:
10-element Array{Int64,1}:
0
0
0
0
0
0
0
0
0
0
In :
I = zeros(T)

Out:
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 :
I = I_0   # Could use OffsetArrays.jl -- enables arbitrary indexing

Out:
1
In :
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 :
I

Out:
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 :
using Plots

┌ Info: Precompiling Plots [91a5bcdd-55d7-5caf-9e0b-520d859cae80]

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

Out:
In :
T = 20

I = zeros(I)
I = 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 

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

## Functions¶

In :
function run_infection(T=20)  # default

I = zeros(T)
I = I_0

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

return I
end

Out:
run_infection (generic function with 2 methods)
In :
methods(run_infection)

Out:
# 2 methods for generic function run_infection:
• run_infection() in Main at In:3
• run_infection(T) in Main at In:3
In :
run_infection(10)

Out:
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 :
I_result = run_infection(10);

In :
I_result

Out:
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 :
I_result = run_infection(20)

plot(I_result, m=:o)

Out:
In :
using Interact


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

## Exponential growth¶

In :
plot(run_infection(1000))

Out:
In :
plot(run_infection(1000), yscale=:log10)

Out:
In :
I_result = run_infection(1000)

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

Out:

$\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 :
p = 0.02
α = 0.01
N = 1000

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

function run_infection(T=20)  # default

I = zeros(T)
I = 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:
run_infection (generic function with 2 methods)
In :
I = run_infection(20)

Out:
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 :
plot(I, m=:o)

Out:
In :
plot(I, m=:o, yscale=:log10)

Out:
In :
I = run_infection(100)
plot(I, m=:o, yscale=:log10)

Out:
In :
p = 0.01
α = 0.1
N = 100

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

Out:

Sigmoid shape -- S-shaped

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

Out:

"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 :
rand()

Out:
0.2767437763441829
In :
rand()

Out:
0.7089008345855983
In :
rand()

Out:
0.3789568897206148

www.random.org

In :
randn()

Out:
-0.33523291258515997
In [ ]:


In :
p = 0.02
α = 0.01
N = 1000

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

I = zeros(T)
I = 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:
run_infection (generic function with 3 methods)
In :
c_average = 1.1

cs = [c_average + 0.1*randn() for i in 1:100]

Out:
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 :
scatter(cs)

Out:
In [ ]: