Remove all objects from workspace.
r remove (list = objects() ) Load add-on packages - deSolve - contains lsoda function - differential equation solver.
library (deSolve)
##
## Attaching package: 'deSolve'
##
## The following object is masked from 'package:graphics':
##
## matplot
Function to compute derivatives of the differential equations.
tiv_model = function (current_timepoint, state_values, parameters)
{
# create state variables (local variables)
T = state_values [1] # target cell
I = state_values [2] # infected cell
V = state_values [3] # virus
with (
as.list (parameters), # variable names within parameters can be used
{
# compute derivatives
dT = lambda - (d * T) - (k *V * T)
dI = (k *V * T) - (delta * I)
dV = (p * I) - (c * V)
# combine results
results = c (dT, dI, dV)
list (results)
}
)
}
Parameters
lambda = 10 #reproduction rate of uninfected cells per cubic millimeter per day
d = 0.01 #death rate of unifected cells per day
k = 0.0000457 #HIV infection rate per cubic millimeter per copies per day
delta = 0.4 #death rate of infected cells per day
p = 38 #virus production rate per copies per cell per day
c = 2.4 # virus clearence rate per day
Disease dynamics parameters.
parameter_list = c (lambda, d, k, delta, p, c)
Initial values for sub-populations.
X = 1000 # target cells
Y = 0 # infected cells
Z = .01 # infectious virus
Compute total cell and viral counts.
N = X + Y + Z
Initial state values for the differential equations.
initial_values = c (T = X, I = Y, V = Z)
Output timepoints.
timepoints = seq (0, 900, by=1)
Simulate the TIV model.
output = lsoda (initial_values, timepoints, tiv_model, parameter_list)
Plot dynamics of target cells.
plot (T ~ time, data = output, type='b', col = 'blue')
Plot dynamics of infected cells.
plot (I ~ time, data = output, type='b', col = 'red')
Plot dynamics of infectious virus.
plot (V ~ time, data = output, type='b', col = 'purple')
Plot dynamics of target cells, infected cells, and infectious virus in the same plot.
# target cells over time
plot (T ~ time, data = output, type='b', ylim = c(0,2000), col = 'blue', ylab = 'T, I, V', main = 'Viral and Immune System Dynamics of HIV Using TIV Model')
text("T", col = 'blue')
## Warning in xy.coords(x, y, recycle = TRUE): NAs introduced by coercion
# remain on same frame
par (new = TRUE)
# infected cells over time
plot (I ~ time, data = output, type='b', ylim = c(0,2000), col = 'red', ylab = '', axes = FALSE)
text("I", col = 'red')
## Warning in xy.coords(x, y, recycle = TRUE): NAs introduced by coercion
# remain on same frame
par (new = TRUE)
# infectious virus over time
plot (V ~ time, data = output, type='b', ylim = c(0,2000), col = 'purple', ylab = '', axes = FALSE)
text("V", col = 'purple')
## Warning in xy.coords(x, y, recycle = TRUE): NAs introduced by coercion
This graph represents a simualtion of the viral and immune system dynamics of HIV with no treatment. At the start of the simulation, the target cell count is 1000, infectious virus (viral load) is 0.1, and no cells have been infected. As we get to about day 50, the viral load has increased along with the infected cells. As a result, the number of target cells have decreased. This trend continues and fluctuates until about day 350 where the viral load, target cells, and infected cells all parallel each other.