2013-12-19 2 views
0

Я пытаюсь решить модель SVEIR (восприимчивая, вакцинированная, подвергнутая воздействию, зараженная и удаленная), используя deSolve. Вспышка начинается на 8-й день (путем импорта индекса в восприимчивой популяции). Для захвата этого я использую событие (путем добавления значения одного (1) к переменной состояния (I) в момент времени Т = 8.Значения разных параметров для разных временных интервалов в deSolve

# Model's parameters 

parms <- c(beta=1.29, 
     betaE=0.25, 
     betaI=1, 
     betaV=0.0, 
     sigma=0.5, 
     gama=0.2, 
     delta=1/365, 
     m=0.000046, 
     r=0.000052, 
     kapa=1.857/10000,    
     alpha=0.00643, 
     thita=1/365, 
     f=0.002)  
dt <- seq(0,50,0.25)  

inits <- c(S=14900, V=0, E=0, I=0, R=0)  
N <- sum(inits) 

eventdat <- data.frame(var = c("I"),time = c(8), 
        value = c(1), method = c("add")) 
eventdat 

#The SVEIR model 

SVEIR <- function(t, x, parms){ 

with(as.list(c(parms,x)),{ 
dS <- - beta*betaE*E*(S/N) - beta*betaI*I*(S/N) - f*S - m*S +delta*R + thita*V + r*N 
dV <- - beta*betaE*betaV*E*(V/N) - beta*betaI*betaV*I*(V/N) - m*V - thita*V + f*S 
dE <- + beta*betaE*E*(S/N) + beta*betaI*I*(S/N) + beta*betaE*betaV*E*(V/N) + beta*betaI*betaV*I*(V/N) - (m + kapa + sigma)*E 
dI <- + sigma*E - (m + alpha + gama)*I 
dR <- kapa*E + gama*I - m*R - delta*R  
der <- c(dS, dV, dE, dI, dR) 
list(der)  
}) 

} 

library(deSolve) 

out <- as.data.frame(lsoda(inits, dt, SVEIR, parms=parms, events = list(data = eventdat))) 

# Plotting the output 

attach(out) 

matplot(x = out[,1], y = out[,-1], type = "l", lwd = 2, 
    lty = "solid", col = c("red", "blue", "black", "green", "darkgreen"), 
    xlab = "time", ylab = "y", main = "SVEIR model") 

legend("bottomright", col = c("red", "blue", "black", "green", "darkgreen"), 
    legend = c("S", "V", "E", "I", "R"), lwd = 2) 

Кроме того, я хочу, чтобы моя модель также захвата . изменения в некоторых параметров Итак, я пытался (пока безуспешно) интегрировать в моей функции а «а» или «за» цикл, который принимает во внимание следующее:

  1. за период времени между 0 - 9 Мне нужно значение параметра betaV быть 0
  2. в течение периода времени между 10 - 50 мне нужно значение betaV параметр, чтобы быть 0,002

Я пытался использовать событие, но R дает мне ошибку (я предполагаю, что я могу использовать событие только для переменных, а не для параметров).

Любая идея, как можно справиться с этим?

Спасибо большое,

Том

PS: (.. Samsuzzoha и др 2012) Модель основана на работе.

ответ

0

Ваш основной вопрос заключается в том, как указать два разных значения betaV в зависимости от времени. Разве вы не можете сделать это в функции, как:

#The SVEIR model 
SVEIR <- function(t, x, parms){ 
    with(as.list(c(parms,x)),{ 
    betaV <- ifelse(t<10,betaV,0.002) # adjust betaV based on value of t 
    dS <- - beta*betaE*E*(S/N) - beta*betaI*I*(S/N) - f*S - m*S +delta*R + thita*V + r*N 
    dV <- - beta*betaE*betaV*E*(V/N) - beta*betaI*betaV*I*(V/N) - m*V - thita*V + f*S 
    dE <- + beta*betaE*E*(S/N) + beta*betaI*I*(S/N) + beta*betaE*betaV*E*(V/N) + beta*betaI*betaV*I*(V/N) - (m + kapa + sigma)*E 
    dI <- + sigma*E - (m + alpha + gama)*I 
    dR <- kapa*E + gama*I - m*R - delta*R  
    der <- c(dS, dV, dE, dI, dR) 
    list(der)  
    }) 

Обратите внимание, что ваш вопрос на самом деле не указать значение для betaV на 9 < t < 10, поэтому я принял обрезание на 10

Когда я запускаю это с betaV = 0.002 (t>10), нет различий в выходе. Если я установил betaV в 1 или 10 для t > 10, V(t) подавлен для больших t и S, E, I and R сдвинут в более низкое время. Правильно ли это звучит?

+0

Уважаемый jlhoward, Теперь он отлично работает. Вы также правы (относительно отсечки на 10). Большое спасибо за вашу ценную помощь. С уважением, Tom – Tom

+0

Добро пожаловать. Рад помочь. Если ответ был полезен, подумайте о «принятии» его ([здесь приведены рекомендации SO] (http://stackoverflow.com/help/someone-answers)). – jlhoward

+0

Сделано !! Еще раз спасибо.... – Tom

Смежные вопросы