Hello! I am having an issue trying to solve differential equations for my model. I had a previous model working but added one component and removed two others and now I cannot get this code to work. I have checked for errors multiple times and I cannot find my issue. Thank you so much.
library(deSolve)
library(plotly)
library(tidyverse)
CCHFModelNew = function(t,x,params)
{
get SIR values
SH = x[1]
EH = x[2]
IH = x[3]
RH = x[4]
ST = x[5]
IT = x[6]
SC = x[7]
IC = x[8]
RC = x[9]
get params
Beta values
betaHH = params["betaHH"]
betaTH = params["betaTH"]
betaCH = params["betaCH"]
betaTC = params["betaTC"]
betaCT = params["betaCT"]
betaTTV = params["betaTTV"] # vertical transmission
betaTTH = params["betaTTH"]
Gamma value
gamma = params["gamma"]
death rates
muH = params["muH"]
muT = params["muT"]
muC = params["muC"]
birth rates
piH = params["piH"]
piT = params["piT"]
piC = params["piC"]
incubation
deltaH = params["deltaH"]
recovery rate
alpha1 = params["alpha1"]
alpha2 = params["alpha2"]
total population
NH = params["NH"] #(SH + IH + EH + RH) + (piH * SH) - (muH * SH)
NT = params["NT"] #(ST + IT) + (piT * ST) - (muT * ST)
NC = params["NC"] #(SC + IC + RC) + (piC * SC) - (muH * SC)
#computations
dSHdt <- piH * NH - ((betaHH * SH * IH) / NH) - ((betaCH * SH * IC) / NH) - ((betaTH * SH * IT) / NH) - (muH * SH)
dEHdt <- (betaHH * SH * IH)/NH + (betaCH * SH * IC)/NH + (betaTH * SH * IT)/NH - ((deltaH + muH)EH)
dIHdt <- (deltaH * EH) - ((alpha1 + gamma + muH) IH)
dRHdt <- (alpha1 * IH) - (muH * RH)
dSTdt <- piT * NT - (betaCT * ST * IC)/NT - (betaTTV * ST * IT)/NT - (betaTTH * ST * IT)/NT - (muT * ST)
dITdt <- (betaCT * ST * IC)/NT + (betaTTV * ST * IT)/NT + (betaTTH * ST * IT)/NT - (muT * IT)
dSCdt <- piC * NC - (betaTC * SC * IT)/NC - (muC * SC)
dICdt <- (betaTC * SC * IT)/NC - ((alpha2 + muC) * IC)
dRCdt <- (alpha2 * IC) - (muC * RC)
return results
list(c(dSHdt, dEHdt, dIHdt, dRHdt, dSTdt, dITdt, dSCdt, dICdt, dRCdt))
}
params = c(betaHH = .096,
betaTH = .157,
betaCH = .067,
betaTC = (1/365), # One tick attaches to one carrier per year
betaCT = 99/365, # One cattle infects 99 ticks per year (assuming 100 ticks on cattle)
betaTTV = ((1/(365 * 2)) * .04) * 280, # ticks give birth once in a lifetime to 7000 ticks with a 4% chance of vertical transmission
betaTTH = 1/365, # check horizontal transmission
gamma = 1/10, # death occurs 7-9th day after onset of illness plus 2 day incubation
muH = (1/(365 * 79)),
muT = (1/(365* 2)),
muC = (1/(8 * 365)), #sheep/deer live 6-11 years
piH = 1.25/(79 * 365), # one couple produces 2.5 children in a lifetime, so one mother produces 1.25
piT = 1 /365, #280/ ( 365 * 2), # 4% of eggs survive
piC = 7/(8 * 365), # sheep produce 7 babies in their life
deltaH = 1/2,
alpha1 = 1/17, # recovery after 15 days
aplha2 = 1/7,
NH = 1100,
NT = 1100,
NC = 1100)
time to start solution
time = seq(from = 0, to = 365, by = 1)
#initialize initial conditions
initialX = c(SH = 1000, EH = 2, IH = 1, RH = 0, ST = 1000, IT = 4, SC = 1000, IC = 1, RC = 0)
CCHFData = ode(y = initialX, times = time, func = CCHFModelNew, parms = params)%>%
as.data.frame()
graph data
data4 <- pivot_longer(data=CCHFData,
cols=-time,
names_to = "initialX",
values_to="value")
myPlot2 <- ggplot(data4, aes(x=time, y= value,color=initialX)) + geom_line()
ggplotly(p = myPlot2, tooltip = "all", dynamicTicks = TRUE, originalData = TRUE)