|
| 1 | +--- |
| 2 | +interact_link: notebooks/ross2010/julia.ipynb |
| 3 | +title: 'Julia' |
| 4 | +permalink: 'chapters/ross2010/julia' |
| 5 | +previouschapter: |
| 6 | + url: chapters/ross2010/intro |
| 7 | + title: 'SIRS dynamics in a large population of households' |
| 8 | +nextchapter: |
| 9 | + url: chapters/network |
| 10 | + title: 'Network models' |
| 11 | +redirect_from: |
| 12 | + - 'chapters/ross2010/julia' |
| 13 | +--- |
| 14 | + |
| 15 | +## Disease dynamics in a population of households in Julia |
| 16 | + |
| 17 | +*Author*: Tim Kinyanjui @timothykinyanjui |
| 18 | + |
| 19 | +*Date*: 2018-10-04 |
| 20 | + |
| 21 | + |
| 22 | +{:.input_area} |
| 23 | +```julia |
| 24 | +# Set the required model parameters for the SIRS model with two levels of transmission - Within and between households |
| 25 | +N = 10; # Household size - Change to 10 for final analysis |
| 26 | +betaHH = 6; # Within household transmission parameter |
| 27 | +betaG = 1; # Population wide transmission |
| 28 | +gamma = 1; # Rate of recovery from infection |
| 29 | +tau = 1; # Rate of loss of protection |
| 30 | +params = [betaHH,gamma,tau,betaG,N]; # Put all the parameters together |
| 31 | +time = (0.0, 30.0) # Simulation time - note it defined as a float |
| 32 | +dim = dim = 0.5*(N+1)*(N+2); # Number of possible configurations - works for three epidemiological classes |
| 33 | +y0 = vec(zeros(1,dim)); # Initial condition vector |
| 34 | +y0[end-1] = 0.00000001; |
| 35 | +y0[end] = 0.99999999; |
| 36 | +``` |
| 37 | + |
| 38 | + |
| 39 | +{:.input_area} |
| 40 | +```julia |
| 41 | +function hhTransitions(N,dim) |
| 42 | + # Function to generate transition matrices for household model |
| 43 | + # Input: N is the household size |
| 44 | + |
| 45 | + # Initialize things |
| 46 | + Qinf = zeros(dim,dim); |
| 47 | + Qrec = zeros(dim,dim); |
| 48 | + Qext = zeros(dim,dim); |
| 49 | + Qwane = zeros(dim,dim); |
| 50 | + dataI = Array{Int64}(zeros(dim,3)) |
| 51 | + m = 0; |
| 52 | + I = Array{Int64}(zeros(N+1,N+1)) |
| 53 | + |
| 54 | + # To help remember where to store the variables |
| 55 | + for ss = 0:N |
| 56 | + for ii = 0:(N-ss) |
| 57 | + m = m + 1; |
| 58 | + I[ss+1,ii+1] = m |
| 59 | + end |
| 60 | + end |
| 61 | + |
| 62 | + # Describe the epidemiological transitions |
| 63 | + |
| 64 | + # Counter for susceptibles |
| 65 | + for ss = 0:N |
| 66 | + # Counter for infecteds |
| 67 | + for ii = 0:(N-ss) |
| 68 | + # If susceptibles and infecteds are more than 1, then infection within the household can occur |
| 69 | + if (ss > 0 && ii > 0) |
| 70 | + Qinf[I[ss+1,ii+1],I[ss,ii+2]] = ii*ss/(N-1); |
| 71 | + end |
| 72 | + |
| 73 | + # If infecteds are more than 1, recovery can occur |
| 74 | + if ii > 0 |
| 75 | + # Rate of recovery |
| 76 | + Qrec[I[ss+1,ii+1],I[ss+1,ii]] = ii; |
| 77 | + end |
| 78 | + |
| 79 | + # For external infection - just keep track of susceptibles |
| 80 | + if ss > 0 |
| 81 | + # Rate of within household infection |
| 82 | + Qext[I[ss+1,ii+1],I[ss,ii+2]] = ss; |
| 83 | + end |
| 84 | + |
| 85 | + # Loss of protection hence becoming susceptible again. Possible if N-ss-ii = rr > 0 |
| 86 | + if (N-ss-ii) > 0 |
| 87 | + # Rate of loss of protection |
| 88 | + Qwane[I[ss+1,ii+1],I[ss+2,ii+1]] = N-ss-ii; |
| 89 | + end |
| 90 | + |
| 91 | + # Store the relevant indices to help identify the household configurations |
| 92 | + dataI[I[ss+1,ii+1],:] = [ss, ii, N-ss-ii]; |
| 93 | + end |
| 94 | + end |
| 95 | + |
| 96 | + Qinf = Qinf - diagm(vec(sum(Qinf,2)),0); |
| 97 | + Qrec = Qrec - diagm(vec(sum(Qrec,2)),0); |
| 98 | + Qext = Qext - diagm(vec(sum(Qext,2)),0); |
| 99 | + Qwane = Qwane - diagm(vec(sum(Qwane,2)),0); |
| 100 | + |
| 101 | + # Return |
| 102 | + return Qinf, Qrec, Qext, Qwane, dataI |
| 103 | +end |
| 104 | +Qinf, Qrec, Qext, Qwane, dataI = hhTransitions(N,dim); |
| 105 | +``` |
| 106 | + |
| 107 | + |
| 108 | +{:.input_area} |
| 109 | +```julia |
| 110 | +using DifferentialEquations |
| 111 | +function rateSIRS(dy_dt,y0,params,time) |
| 112 | + |
| 113 | + # Extract the parameters |
| 114 | + betaHH = params[1]; |
| 115 | + gamma = params[2]; |
| 116 | + tau = params[3]; |
| 117 | + betaG = params[4]; |
| 118 | + N = params[5]; |
| 119 | + |
| 120 | + # Generate the transition matrices |
| 121 | + Qinf, Qrec, Qext, Qwane, HHconfig = hhTransitions(N,dim); |
| 122 | + |
| 123 | + # Combine within and external transitions |
| 124 | + Q = betaHH*Qinf + gamma*Qrec + tau*Qwane + (betaG*((HHconfig[:,2]'*y0)/N)*Qext); |
| 125 | + |
| 126 | + # Generate the rates of change for ODE solver |
| 127 | + dy_dt .= (y0'*Q)'; |
| 128 | + |
| 129 | + # Alternatively this works too |
| 130 | + #=for i=1:length(y0) |
| 131 | + dy_dt[i] = y0'*Q[:,i] |
| 132 | + end=# |
| 133 | + |
| 134 | +end |
| 135 | + |
| 136 | +# Define the ODE problem |
| 137 | +prob = ODEProblem(rateSIRS,y0,time,params); |
| 138 | + |
| 139 | +# Solve |
| 140 | +sol = solve(prob); |
| 141 | +``` |
| 142 | +
|
| 143 | +
|
| 144 | +{:.input_area} |
| 145 | +```julia |
| 146 | +using Plots |
| 147 | +``` |
| 148 | +
|
| 149 | +
|
| 150 | +{:.input_area} |
| 151 | +```julia |
| 152 | +Iconfig = dataI[:,2]; |
| 153 | +infProp = zeros(length(sol.t),1); |
| 154 | + |
| 155 | +# Prepare the plots |
| 156 | +for i = 1:length(sol.t) |
| 157 | + infProp[i,1] = sol[:,i]'*Iconfig/N; |
| 158 | +end |
| 159 | +``` |
| 160 | +
|
| 161 | +
|
| 162 | +{:.input_area} |
| 163 | +```julia |
| 164 | +# Total infectious in the population |
| 165 | +plot(sol.t,infProp,color="blue",xlabel="Time",ylabel="Proportion infectious",label=["Mean infection"],ylims=[0, 1]) |
| 166 | +``` |
| 167 | +
|
| 168 | +
|
| 169 | +
|
| 170 | +
|
| 171 | + |
| 172 | +
|
| 173 | +
|
| 174 | +
|
| 175 | +
|
| 176 | +{:.input_area} |
| 177 | +```julia |
| 178 | +# Household profile at endemic prevalence |
| 179 | +# Prepare the plots |
| 180 | +using PyPlot |
| 181 | +yprop = zeros(N+1,length(sol.t)) |
| 182 | +for j = 1:length(sol.t) |
| 183 | + for i = 1:N+1 |
| 184 | + index = find(Iconfig.==i-1); |
| 185 | + yprop[i,j] = sum(sol[index,j]) |
| 186 | + end |
| 187 | +end |
| 188 | +step(-0.5:1:10.5,[yprop[:,length(sol.t)];0]) |
| 189 | +``` |
| 190 | +
|
| 191 | +
|
| 192 | + |
| 193 | +
|
| 194 | +
|
| 195 | +
|
| 196 | +
|
| 197 | +
|
| 198 | +{:.output_data_text} |
| 199 | +``` |
| 200 | +1-element Array{PyCall.PyObject,1}: |
| 201 | + PyObject <matplotlib.lines.Line2D object at 0x7fe91ef1c210> |
| 202 | +``` |
| 203 | +
|
| 204 | +
|
| 205 | +
|
| 206 | +
|
| 207 | +{:.input_area} |
| 208 | +```julia |
| 209 | +# Household profile at peak prevalence |
| 210 | +# Prepare the plots |
| 211 | +x = find(infProp.==maximum(infProp)) |
| 212 | +step(-0.5:1:10.5,[yprop[:,x];0]); |
| 213 | +``` |
| 214 | +
|
| 215 | +
|
| 216 | + |
| 217 | +
|
0 commit comments