VIRAL AND IMMUNE SYSTEM DYNAMICS OF HIV

Remove all objects from workspace.

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.

T_I_V = function (current_timepoint, state_values, parameters)
{
  # create state variables (Local variables)
  T = state_values [1]        # Uninfcted CD4 T cells (Targets)
  I = state_values [2]        # Infected CD4 T cells
  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) - (sigma * I)
    dV = (p * I) - (c * V)
    
    # combine results
    results = c (dT, dI, dV)
    list (results)
    }
  )
}

Parameters

lambda = 10                   # Reproduction rate of uninfected cells
d = 0.01                      # Death rate of uninfected cells
k = 0.0000457                 # HIV infection rate
sigma = 0.4                   # Death rate of infected cells
p = 38                        # Virus production rate
c = 2.4                       # Virus clearance rate
eRT = 0.4                     # Reverse transcriptase inhibitors effecacy
ePI =  0.4                    # Protease inhibitors efficacy

Disease dynamics parameters

parameter_list = c (lambda, d, k, sigma, p, c)

Initial values for sub-populations.

X = 1000                      # Target cells
Y = 0                         # Infected cells
Z = 0.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, 2000, by=1)

Simulate the TIV model

output = lsoda (initial_values, timepoints, T_I_V, parameter_list)

Plot dynamics of target cells sub-population.

plot (T ~ time, data = output, type='b', col = 'blue')

Plot dynamics of Infected cells sub-population.

plot (I ~ time, data = output, type='b', col = 'red')

Plot dynamics of Infectious viruses sub-population.

plot (V ~ time, data = output, type='b', col = 'green')

Plot dynamics of Targets, Infected cells and Infectious viruses sub-populations 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 = 'TIV epidemic') 

# 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) 

# remain on same frame
par (new = TRUE)  

# infectious viruses over time
plot (V ~ time, data = output, type='b', ylim = c(0,2000), col = 'green', ylab = '', axes = FALSE)

COMMENT:This graph represents a simualtion of the viral and immune system dynamics of HIV when no treatment is applied. The simulation has been run over 2000 days. At the start of the simulation, the target cell count is 1000, infectious virus (viral load) is 0.1, and no CD4 T cells have been infected. We can observe a quick drop in targets cells as the viral load spikes and more cells get infected.

TREATMENT IMPACT

Function to compute derivatives of the differential equations.

T_I_VI_VNI_model = function (current_timepoint, state_values, parameters)
{
  # create state variables (local variables)
  T = state_values [1]        # target cell
  I = state_values [2]        # infected cell
  VI = state_values [3]       # infectious virus 
  VNI = state_values [4]      # noninfectious virus 
  
  
  with ( 
    as.list (parameters),     # variable names within parameters can be used 
         {
           # compute derivatives
           dT = lambda - (d * T) - ((1-eRT) * (k * VI * T))
           dI = ((1-eRT) * (k * VI * T)) - (sigma * I) 
           dVI = ((1-ePI) * (p * I)) - (c * VI)
           dVNI = (ePI * p * I) - (c * VNI)
           
           
           # combine results
           results = c (dT, dI, dVI, dVNI)
           list (results)
         }
    )
}

Parameters

lambda = 10                   # Reproduction rate of uninfected cells
d = 0.01                      # Death rate of uninfected cells
k = 0.0000457                 # HIV infection rate
sigma = 0.4                   # Death rate of infected cells
p = 38                        # Virus production rate
c = 2.4                       # Virus clearance rate
eRT = 0.4                     # Reverse transcriptase inhibitors effecacy
ePI =  0.4                    # Protease inhibitors efficacy

Disease dynamics parameters

parameter_list = c (lambda, d, k, sigma, p, c, eRT, ePI)

Initial values for sub-populations.

W = 0           #noninfectious virus
X = 1000        # target cells
Y = 0           # infected cells
Z = .01         # infectious virus

Compute total cell and viral counts

N = W + X + Y + Z

Initial state values for the differential equations.

initial_values = c (VNI = W, T = X, I = Y, VI = Z)

Output timepoints.

timepoints = seq (0, 1000, by=1)

Simulate the T_I_VI_VNI model.

output = lsoda (initial_values, timepoints, T_I_VI_VNI_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 (VI ~ time, data = output, type='b', col = 'green')

Plot dynamics of noninfectious virus.

plot (VNI ~ time, data = output, type='b', col = 'purple')

Plot dynamics of Targets, Infected cells, Infectious viruses and noninfectious virus sub-populations in the same plot.

# target cells over time
plot (T ~ time, data = output, type='b', ylim = c(0,2000), col = 'blue', ylab = 'T, I, VI, VNI', main = 'T_I_VI_VNI epidemic') 
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 viruses over time
plot (VI ~ time, data = output, type='b', ylim = c(0,2000), col = 'green', ylab = '', axes = FALSE)
text("VI", col = 'green')
## Warning in xy.coords(x, y, recycle = TRUE): NAs introduced by coercion
# remain on same frame
par (new = TRUE)

# noninfectious viruses over time
plot (VNI ~ time, data = output, type='b', ylim = c(0,2000), col = 'purple', ylab = '', axes = FALSE)
text("VNI", col = 'purple')
## Warning in xy.coords(x, y, recycle = TRUE): NAs introduced by coercion

COMMENT: This graph represents the same simualtion when treatment is applied. The simulation has been run over 1000 days. We can observe a faster drop in targets cells as the viral load spikes and quickly drops too. The number of infectious virus also drops as most of the virus load become non infectious.