# EVD transmission model from Khan et al 2015 using Python/PyGOM

Author: Thomas Finnie @twomagpi

Date: 2018-10-03

# housekeeping
import numpy
from pygom import DeterministicOde, TransitionType, SimulateOde, Transition

# first set up the states
states = ['Susceptible_lr',
'Susceptible_hr',
'Exposed',
'Infected',
'Hospitalized',
'Recovered'
]


Table 1

Parameter Description Value
$\phi_H$ Modification parameter for infection rate of high-risk susceptible individuals 1.2-2
$\sigma_I$ Disease-induced death rate of infected individuals 0.10
$\sigma_H$ Disease-induced death rate of hospitalized individuals 0.5
$\theta_I$ Recovery rate of infected individuals 0.1
$\theta_H$ Recovery rate of hospitalized individuals 0.2
$\alpha$ Rate at which latent individuals become infectious 0.1
$\tau$ Hospitalization rate for infected individuals 0.16
$\Pi$ Recruitment rate 1.7
p Fraction of the individuals at high-risk 0.2
$\beta$ Transmission rate of disease Estimated
$1/\mu$ Average life of human 63 years
# now set up the parameters
parameters = ['phi_H',
'sigma_I',
'sigma_H',
'theta_I',
'theta_H',
'alpha',
'tau',
'Rec_rate',
'p',
'beta',
'mu',
'lda',
'eta',
'N'
]

# Set this up as birth, death and transitions
# This could have been done as lifts of the ODE system from the paper but
# being more verbose like this makes us concentrate on the model
transitions = [
Transition(origin='Susceptible_lr',
destination='Exposed',
equation='lda*Susceptible_lr',
transition_type=TransitionType.T),
Transition(origin='Susceptible_hr',
destination='Exposed',
equation='phi_H * lda * Susceptible_hr',
transition_type=TransitionType.T),
Transition(origin='Exposed',
destination='Infected',
equation='alpha * Exposed',
transition_type=TransitionType.T),
Transition(origin='Infected',
destination='Hospitalized',
equation='tau * Infected',
transition_type=TransitionType.T),
Transition(origin='Infected',
destination='Recovered',
equation='theta_I * Infected',
transition_type=TransitionType.T),
Transition(origin='Hospitalized',
destination='Recovered',
equation='theta_H * Infected',
transition_type=TransitionType.T),
]

births_deaths = [
# Births first
Transition(origin='Susceptible_lr',
equation='Rec_rate * (1 - p)',
transition_type=TransitionType.B),
Transition(origin='Susceptible_hr',
equation='Rec_rate * p',
transition_type=TransitionType.B),
# Natural Deaths
Transition(origin='Susceptible_lr',
equation='mu * Susceptible_lr',
transition_type=TransitionType.D),
Transition(origin='Susceptible_hr',
equation='mu * Susceptible_hr',
transition_type=TransitionType.D),
Transition(origin='Exposed',
equation='mu * Exposed',
transition_type=TransitionType.D),
Transition(origin='Infected',
equation='mu * Infected',
transition_type=TransitionType.D),
Transition(origin='Hospitalized',
equation='mu * Hospitalized',
transition_type=TransitionType.D),
Transition(origin='Recovered',
equation='mu * Recovered',
transition_type=TransitionType.D),
#disease deaths
Transition(origin='Infected',
equation='sigma_I * Infected',
transition_type=TransitionType.D),
Transition(origin='Hospitalized',
equation='sigma_H * Hospitalized',
transition_type=TransitionType.D)
]

derived_parameters = [('lda', 'beta * ((Infected + (eta*Hospitalized))/N)'),
('N', 'Susceptible_lr + Susceptible_hr + Exposed + Infected + Hospitalized + Recovered')]

# build the model system
model_system = DeterministicOde(states,
parameters,
transition=transitions,
birth_death=births_deaths,
derived_param=derived_parameters
)
model_system.print_ode()

⎡                                  -(Susceptibleₗᵣ⋅β⋅(Hospitalized⋅η + Infecte
⎢dSusceptible_lr/dt=               ───────────────────────────────────────────
⎢
⎢
⎢                                     -Susceptibleₕᵣ⋅β⋅φ_H⋅(Hospitalized⋅η + I
⎢dSusceptible_hr/dt=                  ────────────────────────────────────────
⎢
⎢
⎢                     -Exposed⋅(α + μ)⋅(Exposed + Hospitalized + Infected + Re
⎢   dExposed/dt=      ────────────────────────────────────────────────────────
⎢
⎢
⎢   dInfected/dt=
⎢
⎢ dHospitalized/dt=
⎢
⎣  dRecovered/dt=

d) + (Recᵣₐₜₑ⋅(p - 1) + Susceptibleₗᵣ⋅μ)⋅(Exposed + Hospitalized + Infected +
──────────────────────────────────────────────────────────────────────────────
Exposed + Hospitalized + Infected + Recovered + Susceptibleₕᵣ + Susceptibleₗᵣ

nfected) + (Recᵣₐₜₑ⋅p - Susceptibleₕᵣ⋅μ)⋅(Exposed + Hospitalized + Infected +
──────────────────────────────────────────────────────────────────────────────
Exposed + Hospitalized + Infected + Recovered + Susceptibleₕᵣ + Susceptibleₗᵣ

covered + Susceptibleₕᵣ + Susceptibleₗᵣ) + Susceptibleₕᵣ⋅β⋅φ_H⋅(Hospitalized⋅η
──────────────────────────────────────────────────────────────────────────────
Exposed + Hospitalized + Infected + Recovered + Susceptibleₕᵣ + Susceptibleₗᵣ

Exposed⋅α - Infected⋅μ - Infected⋅σ_I - Infected⋅τ - Infected⋅θ_I

-Hospitalized⋅μ - Hospitalized⋅σ_H + Infected⋅τ - Infected⋅θ_H

Infected⋅θ_H + Infected⋅θ_I - Recovered⋅μ

Recovered + Susceptibleₕᵣ + Susceptibleₗᵣ))               ⎤
────────────────────────────────────────────              ⎥
⎥
⎥
Recovered + Susceptibleₕᵣ + Susceptibleₗᵣ)                ⎥
──────────────────────────────────────────                ⎥
⎥
⎥
+ Infected) + Susceptibleₗᵣ⋅β⋅(Hospitalized⋅η + Infected)⎥
──────────────────────────────────────────────────────────⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎦


# Create the timeline
t = numpy.linspace(0, 200, 200)
# Set the inital values for the states
model_system.initial_values = ([1000000,
20000,
15,
10,
0,
0],
t[0])

liberia_parameters = {
'phi_H': 1.6,
'sigma_I': 0.1,
'sigma_H': 0.5,
'theta_I': 0.1,
'theta_H': 0.2,
'alpha': 0.1,
'tau': 0.16,
'Rec_rate': 1.7,
'p': 0.2,
'beta': 0.371,
'mu': (1/63)/365,
'eta': 0.7,
}
sierra_leone_parameters = {
'phi_H': 1.6,
'sigma_I': 0.1,
'sigma_H': 0.5,
'theta_I': 0.1,
'theta_H': 0.2,
'alpha': 0.1,
'tau': 0.16,
'Rec_rate': 1.7,
'p': 0.2,
'beta': 0.371,
'mu': (1/63)/365,
'eta': 0.7,
}

model_system.parameters = liberia_parameters
solution,output = model_system.integrate(t[1::], full_output=True)
model_system.plot()

<Figure size 640x480 with 6 Axes>

model_system.parameters = sierra_leone_parameters
solution,output = model_system.integrate(t[1::], full_output=True)
model_system.plot()