Modeling Viral Infections

1. Introduction

Many diseases are caused by viruses.  A viral disease begins with an infection of the body by a small number of viruses.  The viruses attempt to enter individual cells, and "hijack" them into producing new copies of the virus.  An infected cell will die, releasing the new copies of the virus to continue the infection.  If the infection is not stopped, so many cells will be destroyed that the body will die.

The immune system constantly examines the body for "alien" objects, such as viruses.  When it notices an invading virus, it records the pattern of "antigens" on the surface of the virus, and then produces special cells called "antibodies", which are programmed to  destroy all objects that have the same pattern of antigens.  Thus, once the immune system has spotted a single virus, it will try to kill all copies of the virus.

Thus, a typical viral disease usually has just two patterns: either the virus replicates so quickly that the host body is killed, or the immune system is able to kill all the viruses, and the body returns to health.  In either case, the resolution is often quick, a matter of days or weeks.

2.  Modeling a Simple Infection

Let's try to describe the changes that occur in an infection over, say, a single hour.  We can measure the strength of the infection by determining the density of viruses, a number which we will call v(t).   We will also want to measure the density of antibodies produced by the immune system against the virus, which we will call a(t).

v(t)  - number of viruses in the body at time t
  a(t)  - number of antibodies at time t

Our first measurement will occur at time t0  and will give us the values of v(t0) and a(t0). Conditions in the body will cause these values to change with time.  a(t) will rise as the body tries to fight the infection, and 
v(t) will rise as the virus replicates, or fall if the antibodies are able to kill the viruses faster than they reproduce.

What we would like to do is write down rules that summarize this behavior, and allow us to make simple calculations about the course of the infection, without having to take blood samples from a patient every hour.  What we are trying to do now is to model the infection; that is, to try to understand the process of an infection by comparing it to some set of mathematical equations.

First of all, we will assume that the immune system responds to a virus by increasing the number of antibodies in a linear way; that is, the number of antibodies produced is proportional to the number of viruses present. This is a reasonable assumption which we have seen in other models. This results in the following relationship:

      Increase in antibodies in time Δt is k * v(t) * Δt,
  
where k is some number we would have to find experimentally, and Δt is a small interval of time, during which v(t) doesn't change much. For now, we will assume that the antibodies are "immortal", so their number only grows.

Now we need to look at what the virus particles are doing. Each virus particle tries to invade a body cell and force it to make new copies of the virus.  We will suppose that we understand this process completely; that is, we can estimate the chances that a virus will be able to invade a cell at
any particular moment, we know how long the invaded cell will be able to remain alive  making copies, and we know how many copies the cell makes before dying.  If we really know all of that, we can estimate the relative rate r at which new viruses are "born" over a given time interval.  The actual number that are born must be the product of this rate, times the number of viruses, times the amount of time that passes:

     Increase in viruses in time Δt is   r * v(t) * Δt.
  
The only other thing we have ignored so far is that antibodies are busy killing viruses.  Now the antibodies and viruses are wandering around the body, almost at random.  Suppose we know the chances that in time Δt, one particular antibody will run into one particular virus and kill it, and say this chance is p.  Then, the chances that our same antibody will run into any virus in that time must be p*v(t).  And since every antibody has the same chance as the one we picked to run into some virus, the total number of "encounters" is also multiplied by the number of antibodies in the system.  In other words,

    Decrease in viruses in time Δt is    p * v(t) * a(t) * Δt.

Thus, if we knew the values of k, r and p, we could estimate the level of antibodies and viruses in the body  Δt hours  after a measurement at time t by writing:

a(t+Δt) = a(t) + k * v(t) * Δt
  
v(t+Δt) = v(t) + r * v(t) * Δt - p * v(t) * a(t) * Δt

In particular, if we wanted the level of antibodies and viruses in the body 4 hours  after our first  measurement we would simply set t=t0=0 and Δt=4. 

The first equation says that the new number of antibodies is the sum of the old value, plus the newly created antibodies.   The second equation says that the new number of viruses is the sum of the previous number plus the increase in viruses minus the number of viruses destroyed by the antibodies.

Example 1.  Let's look at a simple example of this model.  First, we will assume that some helpful lab technician has figured out the values of the constants k, r and p for a particular virus:

Clear[k,r,p,a,v,delt,t]
k=0.1;
r=0.5;
p=0.5;

Then, we will assume that a patient has just been infected with the virus.   This means that the immune system is not alerted at all, so a(t) will start out as zero; i.e., a(0)=0 where we have set t0=0.  We will give v(t) a small initial value, which is less than 1.  This is because we will assume that v(t) measures a virus density, rather than an actual count of the number of virus particles.  We will use the same convention for a(t), i.e., a(t) will represent the antibody density. For this computation, we will be taking time steps Δt of 12 hours.

a[0]=0.;
v[0]=0.01;
delt=12;

And now, we write down the equations that are supposed to tell us the density of antibodies and viruses 12 hours after time t=0, i.e., when t+Δt=12.

t=0;
a[t+delt] = a[t] + k*v[t]*delt
v[t+delt] = v[t] + r*v[t]*delt - p*v[t]*a[t]*delt

We can repeat this process to find what's happening 24 hours after the initial infection.  Notice that we use the exact same equations as before except that now we know the solution at t=12 and want the solution at t=24.

t=12;
a[t+delt] = a[t] + k*v[t]*delt
v[t+delt] = v[t] + r*v[t]*delt - p*v[t]*a[t]*delt

So at the end of 24 hours the density of the antibodies is 0.096
and the density of the virus has increased to 0.48496.

Problem 1.  

Continue the above calculation for t=36, 48, and so on in intervals of 12 hours. Stop when either the virus density exceeds 20 (which we will assume means that the patient died of the disease) or when the density of the virus becomes zero or negative. (After the density of the virus becomes negative the results are meaningless.)

Print a table of the values of t, a(t) and v(t) using table headings. 

Make a plot of a(t) versus t and v(t) versus t on the same graph. (Use PlotJoined->True.)  Plot only values of t until the stopping criteria is reached.  Use a different color  for each function and label your graph. 

Did the patient survive or did the virus kill him?  After how many hours was the stopping criteria satisfied; i.e., when did the virus density exceed 20 or become negative?

Note: At this point do not define a function for a[t] and v[t] but merely proceed as we did in Example 1. It is possible to define a function for a[t] and v[t] recursively, i.e., in terms a[t-Δt] and v[t-Δt] but the calculations become very slow since Mathematica does not save the values at the previous times but rather continues to recalculate them.  Below we will see how to automate this process.

We now turn to the problem of automating the above procedure. 

Mathematica has a Do command that allows you to make a loop that repeats one or more operations.  The format of the Do command is
 
  Do [ { command1,
         command2,
         ...
         last command } ,
        
        { counter, start-value, end-value, increment} ]
  
For instance, the following Do loop calculates the square of the odd numbers between 1 and 9 and then prints out the number and its square.

Do [ { nsquare = n^2,
        Print ["n=", n, "   n squared= ",nsquare]},
        {n, 1, 9, 2} ]

In this Do loop, the commands to be repeated are grouped inside the curly brackets.  (If there is only one command, you don't need the brackets).  After the commands comes an iterator, similar to what you've seen for the Table command.  This tells Mathematica 

  the variable to be iterated (n in our case)

  how to start the loop (set n to 1 in our case), 
  
  when to stop the loop (n>9 in our case)
    
  how to choose the next value of the iterator (add 2 to the
       previous value in our case).

Example 2. We want to use the Do command to recompute the values of a(t) and v(t) that we got in Example 1.  To be safe we will Clear all  data before starting.

Clear[k,r,p,a,v,t,delt]

k=0.1;
r=0.5;
p=0.5;

delt=12;

a[0]=0.;
v[0]=0.01;

Do [ { a[t+delt] = a[t] + k*v[t]*delt,
        v[t+delt] = v[t] + r*v[t]*delt
                     - p*v[t]*a[t]*delt },
                     {t,0,24,12}]

Note that there is no output from the Do command because we did not ask Mathematica  to print out any results; Mathematica simply performed the calculations.  We could include a Print statement in the Do loop or we could just print out another table since Mathematica  has saved the values.

table2=Table[{t,a[t],v[t]},{t,0,24,12}];

TableForm[table2,TableHeadings->{ None,{"t in hours",
"a(t)","v(t)"}}]

Problem 2.  

Use the Do command to regenerate your table from Problem 1.

Problem 3.  

In the first problems we chose Δt=12 hours.  This was arbitrary and we can choose different intervals to calculate v(t) and a(t). 

Use the Do command to generate a(t) and v(t) at 6 hour intervals.  Make a table of your results and plot them as in Problem 1.   After how many hours was the stopping criteria satisfied; i.e., when did the virus density exceed 20 or become negative?

Use the Do command to generate a(t) and v(t) at 3 hour intervals.  Make a table of your results and plot them as in Problem 1.   After how many hours was the stopping criteria satisfied; i.e., when did the virus density exceed 20 or become negative?
  
Discuss any differences in your results from this problem and Problem 2.

We now turn to the problem of investigating the effect the parameter  r   has on the outcome of the viral infection.  Remember that we got these parameters from experiment and they may include significant experimental errors.  What we would like to do now is investigate what happens in the infection when those numbers change. The following problem investigates this issue.

Problem 4.  

We would like to see whether the virus density can reach the value 20 or greater (killing the patient) before the immune system has a chance to destroy the virus.

Copy your input from Problem 3 using Δt=6 hours in all calculations for this problem.  

Change your Do loop so that it always computes values for 48 hours, which should be enough time to see whether the infection or immune system "wins".  Note that the results are meaningless once any value becomes negative, or any value becomes really large.

The value r determines how fast the virus reproduces.  We started with r=0.5, and the resulting infection died out. What is the smallest value of r (increment r in steps of Δr=.1) that will cause the infection to kill the  patient?  Experiment until you find such an r, and then print a table of results for that value and plot your data as you did in Problem 1.