forked from ICI3D/RTutorials
-
Notifications
You must be signed in to change notification settings - Fork 0
/
particCodingDynMod2016.R
124 lines (107 loc) · 4.32 KB
/
particCodingDynMod2016.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
## How does increased susceptiblility to TB
## amongst diabetics affect incidence of active
## TB?
require(deSolve)
## Big population
## SLAT
slatDB <- function(t,y,parms, browse=F){
with(c(as.list(y),parms),{
if(browse) browser()
N <- S+L+A+T+SDB+LDB+ADB+TDB
TotBirthRate <- birthRate*N + tbDeathRate*A + tbDeathRate*ADB
foi <- transmissionCoeff*(A+ADB)/N ## force of infection
## non-diabetics
dSdt <- -foi*S + clearanceRate*T + TotBirthRate*(1-diabPrev) -
naturalDeathRate*S
dLdt <- foi*S - activationRate*L -
naturalDeathRate*L
dAdt <- activationRate*L - treatmentRate*A - tbDeathRate*A -
naturalDeathRate*A
dTdt <- treatmentRate*A - clearanceRate*T - naturalDeathRate*T
dSdt+dLdt+dAdt+dTdt
## diabetics
dSdtDB <- -foi*SDB + clearanceRate*TDB + TotBirthRate*diabPrev -
naturalDeathRate*SDB
dLdtDB <- foi*SDB - activationRate*idrActivationDB*LDB -
naturalDeathRate*LDB
dAdtDB <- activationRate*idrActivationDB*LDB - treatmentRate*ADB -
tbDeathRate*ADB - naturalDeathRate*ADB
dTdtDB <- treatmentRate*ADB - clearanceRate*TDB - naturalDeathRate*TDB
## dSdt+dLdt+dAdt+dTdt+dSdtDB+dLdtDB+dAdtDB+dTdtDB
list(c(dSdt,dLdt,dAdt,dTdt,dSdtDB,dLdtDB,dAdtDB,dTdtDB))
})
}
params <- c(
transmissionCoeff = 10 ## [per year]
, activationRate = 1/2
, treatmentRate = 1/0.5
, clearanceRate = 1/0.5
, tbDeathRate = 1/3
, birthRate = 1/60
, naturalDeathRate = 1/60
, idrActivationDB = 3 ## [unitless]
, diabPrev = 0.12 ## proportion
, N0 = 52*10^6) ## [people]
params
## activationRate*idrActivationDB
time <- 0
initCond <- params['N0'] * c(
S = 2/3
, L = 1/3
, A = 0
, T = 0)
initCondDB <- c(initCond*(1-params['diabPrev']), initCond*params['diabPrev'])
names(initCondDB)[5:8] <- paste0(names(initCondDB)[5:8], 'DB')
slatDB(t=time,y=initCondDB,parms=params, browse=F)
args(slatDB)
timeseq <- seq(1950, 1970, by = .02)
params['transmissionCoeff'] <- 7
simdat <- data.frame(lsoda(
y = initCondDB, # Initial conditions for population
times = timeseq, # Timepoints for evaluation
func = slatDB, # Function to evaluate
parms = params# Vector of parameters
))
simdat <- within(simdat, { N <- S+L+A+T+SDB+LDB+ADB+TDB})
simdat <- within(simdat, { PA <- (A+ADB)/N})
head(simdat); tail(simdat)
par('ps'=16, lwd = 5, mfrow=c(1,2))
## non-diab
plot(0,0, type = 'n', xlim = range(timeseq), ylim = c(0, params['N0']),
xlab = 'time', ylab = 'people', bty='n')
with(simdat, lines(time, S, col ='black'))
with(simdat, lines(time, L, col ='orange'))
with(simdat, lines(time, A, col ='red'))
with(simdat, lines(time, T, col ='blue'))
legend('topright', leg = c('S','L','A','T'), col = c('black','orange','red','blue'), lwd = 5, cex = 1.5)
## diab
plot(0,0, type = 'n', xlim = range(timeseq), ylim = c(0, params['N0']),
xlab = 'time', ylab = 'people', bty='n')
with(simdat, lines(time, SDB, col ='black'))
with(simdat, lines(time, LDB, col ='orange'))
with(simdat, lines(time, ADB, col ='red'))
with(simdat, lines(time, TDB, col ='blue'))
legend('topright', leg = c('SDB','LDB','ADB','TDB'), col = c('black','orange','red','blue'), lwd = 5, cex = 1.5)
par('ps'=16, lwd = 5, mfrow=c(1,2))
## non-diab
with(simdat, plot(time,PA, type='l'))
## compare PA for differennt levels of diabetes prevalence
diabPrevVec <- seq(0, .4, by = .05)
eqPA <- rep(NA, length(diabPrevVec))
for(jj in 1:length(diabPrevVec)) {
params['transmissionCoeff'] <- 7
params['diabPrev'] <- diabPrevVec[jj]
simdat <- data.frame(lsoda(
y = initCondDB, # Initial conditions for population
times = timeseq, # Timepoints for evaluation
func = slatDB, # Function to evaluate
parms = params# Vector of parameters
))
simdat <- within(simdat, { N <- S+L+A+T+SDB+LDB+ADB+TDB})
simdat <- within(simdat, { PA <- (A+ADB)/N})
eqPA[jj] <- tail(simdat$PA,1)
}
par(cex=1.5, lwd = 3, bty = 'n')
plot(diabPrevVec, eqPA, type='l', xlab = 'diabetes prevalence',
ylab = 'equilibrium prevalence of active TB',
ylim = c(0, max(eqPA)))