# One Host SEIR + One Vector SEI model in Julia

Author: Carl A. B. Pearson
Date: 2 OCT 2018

First, create a system of ODEs conforming to the required signature:

function F(du,u,p,t)
S_H, E_H, I_H, R_H, S_V, E_V, I_V = u

# host dynamics
host_infection = (p.β*S_H*I_V)/p.N_H
host_mortality = p.μ_H .* u[1:4] # include S_H, so easier to remove mortality
host_births = sum(host_mortality)
host_progression = p.σ_H*E_H
recovery = p.λ*I_H

du[1] = -host_infection + host_births
du[2] = host_infection - host_progression
du[3] = host_progression - recovery
du[4] = recovery
du[1:4] -= host_mortality

# vector dynamics
vec_infection = (p.β*S_V*I_H)/p.N_H
vec_mortality = p.μ_V .* u[5:7] # include S_V, so easier to remove mortality
vec_births = sum(vec_mortality)
vec_progression = p.σ_V*E_V

du[5] = -vec_infection + vec_births
du[6] = vec_infection - vec_progression
du[7] = vec_progression
du[5:7] -= vec_mortality

end
F (generic function with 1 method)

Set initial conditions and dynamic parameters, then apply the ODE solver. Note this will also display the solver run time; this will be slow for the first run, though quite fast if you re-run it (even after changing parameters or initial conditions).

using DifferentialEquations
using IterableTables, DataFrames
using NamedTuples
# nb: in >= Julia v0.7, can eliminate this import
#  and the @NT syntax
u0 = [
S_H=100.0,   E_H=0.0, I_H=1.0, R_H=0.0,
S_V=10000.0, E_V=0.0, I_V=0.0
]
p = @NT(
μ_H=1/365, μ_V=1/30, σ_H=1/3, σ_V=1/7, λ=1/14,
β=0.05, N_H = sum(u0[1:4])
)
tspan = (0.0, 365.0)
prob = ODEProblem(F, u0, tspan, p)
sol = @time solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8,saveat=linspace(0,365,365*10+1))

df = DataFrame(sol)
rename!(df, Dict(:timestamp => :t,
:value1 => :S_H, :value2 => :E_H, :value3 => :I_H, :value4 => :R_H,
:value5 => :S_V, :value6 => :E_V, :value7 => :I_V
))
mlt = melt(df,:t) # convert results into long format for plotting
mlt[:host] = contains.(string.(mlt[:variable]),"H"); # tag which entries are host vs vector
df
3.965556 seconds (3.46 M allocations: 192.139 MiB, 1.64% gc time)

tS_HE_HI_HR_HS_VE_VI_V
10.0100.00.01.00.010000.00.00.0
20.1100.000268141649855.745401288701111e-60.99261065501874190.0071154579301127969999.507616077060.48888117779078830.0035027451477435017
30.2100.000501778290164.525792070955373e-50.98527657398145150.0141763898076934799999.0205329147870.96558801279815520.013879072414963943
40.3100.000667379846120.000150405677417930470.97799900538715150.0211832090893128549998.5387046862051.43036111437824620.030934199416335188
50.4100.000732379486280.00035106529563958110.97078021351299320.0281363417050963769998.0620848577951.88343774323796450.05447739896513851
60.5100.000665153699670.00067520708974060010.9636234061566550.035036233053943399997.59062570772.32505235061697270.08432194168190127
70.6100.000435002734850.00114897706930118030.95653266570568070.0418833544901830249997.1242778807172.75543707307972640.12028504620147329
80.7100.000012131410860.00179677567863891240.94951288360361430.048678209306902479996.6629899795973.17482218393747840.16218783646455837
90.899.999367630120190.00264133407201210850.94256969750693880.055421338300874099996.2067081865793.58343650852039540.20985530490011298
100.999.998473455944630.003703787975066140.9357094311792070.062113324901110059995.7553759165283.98150780185506030.2631162816177185
111.099.997302414557980.005003746744429860.92893903895349030.068754799744107779995.3089335091264.3692630814722140.32180340940160024
121.199.99582814114040.0065593644357901440.92226604926203830.075346445161792139994.8673179286494.7469289497420290.38575312160914144
131.299.994025083042430.0083874020535659770.91569851607453990.081888998829476359994.430462522145.114731850682020.45480562717867207
141.399.991868480628370.0105032953036620250.90924496644728310.088383257620694179993.9982967648215.47289833934290.528804895835597
151.499.989334350435770.0129212112109171360.90291435788006720.094830080473251469993.570746073815.8216552755644820.6075986506255675
161.599.986399465764010.015654113011603720.89671602957349840.101230391650906729993.1477315962536.1612300444223090.6910383593243803
171.699.983041340532520.018713810878424380.89065966568124720.107585182907812089992.729170072326.4918506953174210.7789792323625186
181.799.979238209205280.0221110248882701030.88475524971742510.113895516189046719992.3149736639976.8137461163463480.8712802196570351
191.899.97496901140760.0258554305034490940.87901303283279470.120162525256171959991.9050498554647.1271461328974640.9678040116372855
201.999.970213371566390.0299557177089795240.87344349302597930.126387417698668329991.4993013310167.4322816278764391.0684170411069502
212.099.96495158289340.034419636331622850.86805730443147210.13257147634351119991.0976258911617.72938462317628351.1729894856613798
222.199.959164588941990.039254044807858510.86286530587254150.138716060377634879990.6999163913628.0186883355745981.2813952730626423
232.299.952833964559620.0444649616967721650.85787846703885160.14482260670477239990.3060606645858.3004272500464411.3935120853669343
242.399.945941900644020.050057602467425910.8531078664912430.150892630397327739989.9159415039748.574837131139981.5092213648849446
252.499.93847118250940.056036433414422740.84856465845536290.156927725620829029989.529436618028.8421550611546391.6284083208239846
262.599.930405174983550.062405206828083330.84426005218002810.162929566008354629989.146418618049.1026194473815441.7509619345776433
272.699.921727803542080.069167002685000870.84020528899026470.168899904782672229988.766755026769.3564700067263691.8767749665126119
282.799.912423534056770.076324275188115520.83641161560847530.174840575146654059988.3903082692259.6039477681509242.00574396262303
292.899.902477358774680.083878880923612760.83289027025383510.180753490047886359988.0169357016689.8452950392002722.1377692591325457
302.999.891874774188370.09183212455953080.82965245924674660.186640642005375459987.64648964299210.0807553680567212.2727549889506884

Now that we have a solution, we want to view what is happening in host vs mosquito population: