Skip to content

Instantly share code, notes, and snippets.

@sathomas

sathomas/README.md

Last active August 29, 2015 14:16
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save sathomas/3519667ff6b96d4d7c30 to your computer and use it in GitHub Desktop.
Save sathomas/3519667ff6b96d4d7c30 to your computer and use it in GitHub Desktop.
Download ZIP
Mathematical Models of HIV, Part 1
Raw
README.md

Mathematical models of disease often help us understand the epidemiology of disease—how a disease epidemic spreads within a population and how measures such as vaccination can control that spread. This post, however, considers the virology and immunology of HIV rather than its epidemiology. The discussion that follows focuses on the actions of HIV within a single human host.

Raw
README2.md

Mathematical models of disease often help us understand the epidemiology of disease—how a disease epidemic spreads within a population and how measures such as vaccination can control that spread. This post, however, considers the virology and immunology of HIV rather than its epidemiology. The discussion that follows focuses on the actions of HIV within a single human host.

Before Infection

As a starting point we can model the dynamics within a human host before any infection occurs with a particular focus on a type of white blood cell known as a CD4+ T cell. This cell is a favorite target of the virus. (The fact that HIV targets CD4+ T cells is one of the reasons that AIDS is such a devastating disease; CD4+ T cells play a critical role in the immune system’s response to any pathogen.) In this model we assume that the body produces CD4+ T cells at a (roughly) constant rate λ and that the cells have an average lifetime of 1/d. Those assumptions result in a simple model for the dynamics of healthy CD4+ T cells. We can represent the concentration of these cells in the blood using the variable x. (We’re using a dot above a variable name (e.g. ẋ) to represent the derivative of that variable with respect to time (dx/dt).)

$ dot x = lambda - d * x $

This model is a linear differential equation with a straightforward solution for x as a function of time.

$ x(t) = lambda/d + Ce ^ (-d * t) $

The exponential term decays to 0, so that the model shows a (relatively) stable concentration of CD4+ T cells. With a couple of additional data points, we can estimate values for the parameters λ and d.

  • [Nowak] notes that a healthy human adult blood contains CD4+ T cells in a concentration of about 1000/µl.
  • Using clinical data, [Stafford] derives an estimate for d of about 0.01/day.

Initial Infection

When HIV invades the body, it changes the dynamics of CD4+ T cells. We’ll model those changes by introducing two new parameters, y, and v. The quantity y represents CD4+ T cells that have been infected by the virus; we’ll reserve x for healthy CD4+ T cells that have not been infected. The quantity v represents the number of free virus particles (virions) in the body.

In the initial model, x increases when the body produces new cells and decreases when those cells reach the end of their lives. With the introduction of HIV, there is another factor that reduces healthy cells: those cells may become infected with the virus. The rate at which that happens depends on three factors:

  1. The amount of healthy cells that are available for the virus to infect.
  2. The amount of virions circulating in the body that can infect a cell.
  3. The effectiveness of the virus in evading the body’s immune system and infecting a healthy cell.

We represent that third element by the parameter β and update the model’s accounting of healthy CD4+ T cells.

$ dot x = lambda - d x - beta x v $

The βxv factor represents both a decrease in the amount of healthy cells and in increase in the amount of infected cells. Infected cells also have an average lifetime, and it’s plausible to assume that the lifetime of infected cells differs from that of healthy cells. We’ll use a (instead of d) to model the death of infected cells, so that the equation for infected cells becomes:

$ dot y = beta x v - a y $

The final component of our model considers the dynamics of free virus particles. When the virus infects a cell, it hijacks the cell’s reproductive mechanisms for its own reproduction, turning the cell into a virus factory. We define the parameter k to represent the rate at which infected cells produce new virus particles. For the average lifetime of a free virus particle we’ll use the parameter 1/u. The third equation in the model is:

$ dot v = k y - u v $

Our model then becomes a system of three differential equations.

$ [ [ dot x ] , [ dot y ], [ dot v ] ] = [ [ lambda - d x - beta x v ], [ beta x v - a y ], [ k y - u v ] ] $

Unfortunately, the βxv factor in the first two equations is the product of system variables, not the sum. As a result, the system of equations is no longer linear and does not have a known analytical solution. We can still use the model, however, to help understand the course of an HIV infection.

Note: To be clear, this model is a significant simplification of the actual virus dynamics. As we’ll see, though, it can still provide useful insights.

Equilibrium Points

Even though we can’t write a set of equations that describe how the system (x, y, and v) evolves in time, we can find the system’s equilibrium points and determine their stability. Equilibrium points occur when all of the derivatives are zero. In this state the system does not change. (That’s true regardless of the point’s stability.)

Stability determines how the system behaves when it is close to (but not exactly at) an equilibrium point. If the system reaches a state close to a stable equilibrium point, then it will, on its own, evolve towards that equilibrium point over time. If, on the other hand, the system reaches a state close to an unstable equilibrium point, the system will evolve away from that point. Stable equilibrium points are like magnets that attract, while unstable equilibrium points repel.

Finding the equilibrium points for our model requires only algebra. Specifically, we have to solve the system of equations:

$ [ [ 0 ] , [ 0 ], [ 0 ] ] = [ [ lambda - d x - beta x v ], [ beta x v - a y ], [ k y - u v ] ] $

Without boring you with the details, it turns out that our model has two equilibrium points. We’ll consider each separately.

Disease-Free Equilibrium

One equilibrium is not surprising. It is the point (x, y, v) = (λ/d, 0, 0). This point corresponds to the system with no infected cells and no free virus particles. It’s the same as if the virus had never entered the system at all. We call this the disease-free equilibrium. If the system ever reaches this point exactly, then it remains there indefinitely. (Unfortunately, this is true only in our greatly simplified model; in fact, HIV can persist in the body even when there are no free virus particles nor infected CD4+ T cells.)

An important question is whether or not this equilibrium point is stable. If the disease-free equilibrium is stable, then a successful HIV treatment need only move the system “close enough” to the equilibrium. Once the patient is “close enough” to disease-free, the system itself will take over and tend toward the disease-free equilibrium on its own. If the equilibrium is not stable, however, then “close enough” won’t be good enough. No matter how close a treatment forces the system towards the disease-free state, the system’s internal dynamics will resist further movement towards that state.

The stability of the disease-free equilibrium can be determined using linear stability analysis. We’ll skip the details (TLDR: examine the eigenvalues of the system’s Jacobean) and consider the results. It turns out that the stability depends on the values of the system parameters. That stability can be most conveniently captured by combining those parameters into a single new parameter R0.

$ R_0 = ( beta k lambda ) / ( d a u ) $

If R0 is less than 1, then the disease-free equilibrium is stable. If it is greater than 1, then the equilibrium point is not stable. To see how different values of R0 affect the system behavior, consider the visualization that headlines this article.

The graph clearly shows how dramatic a difference the value of R0 makes. If R0 is less than 1, then the number of virus particles tends to zero, even when the system is initialized with a billion free virus particles. If, however, R0 is greater than 1, the system tends towards a state with millions of virus particles in the body, even when the infection bgeins with just a single virus particle.

If R0 for HIV were less than 1, then the system dynamics could naturally drive the system towards the disease-free equilibrium point, and HIV patients would be cured. Sadly, the lowest estimates (from [Nowak]) for HIV’s R0 are at least 5, suggesting that the disease-free equilibrium point is not stable.

Part 2 of this discussion continues by considering the second equilibrium point.

References

Disclaimer

I’m neither a virologist nor a mathematician, just a geek that likes math. Undoubtedly some of what I’ve said above is wrong. When you notice it, please let me know, and I’ll do my best to fix it.

Raw
index.html
<!DOCTYPE html>
<html>
<head>
<meta charset='utf-8'>
<title>Mathematical Models of HIV, Part 1</title>
<link href='http://fonts.googleapis.com/css?family=Varela'
rel='stylesheet' type='text/css'>
<style>
body { font-family: Varela,sans-serif; color: #444; }
</style>
</head>
<body>
<script src='http://d3js.org/d3.v3.min.js'></script>
<script>
// Define the dimensions of the visualization. The
// width and height dimensons are conventional for
// visualizations displayed on http://jsDataV.is.
var margin = {top: 20, right: 10, bottom: 48, left: 60},
width = 636 - margin.left - margin.right,
height = 476 - margin.top - margin.bottom;
// Create the SVG stage for the visualization and
// define its dimensions. Web browsers don't require
// the full SVG attributes, but adding them makes it
// easier to extract the SVG from a web page and
// insert it into a stand-alone file.
var svg = d3.select("body").append("svg")
.attr("width", width + margin.left + margin.right)
.attr("height", height + margin.top + margin.bottom)
.attr("viewBox", "0 0 " +
(width + margin.left + margin.right) + " " +
(height + margin.top + margin.bottom) )
.attr("version", "1.1")
.attr("xmlns", "http://www.w3.org/2000/svg");
// Within the SVG container, add a group element (<g>)
// that can be transformed via a translation to account
// for the margins. This element is restricted to the
// clipping path.
var g = svg.append("g")
.attr("transform", "translate(" + margin.left +
"," + margin.top + ")");
// Parameters for the simulated system. Explanations
// for these parameters may be found in the associated
// text (README.md).
var d = 0.01, // death rate of healthy cells
a = 1.00, // death rate of infected cells
u = 23.0, // death rate of free virus particles
k = 1000, // rate of virus generation
lambda = 100000, // rate of healthy cell generation
beta = 2e-8; // viral efficiency
// The derived parameters that the simulation actually
// used.
var alpha = a / d, // α = a/d
epsilon = d / u; // ε = d/u
// Although we're going to try different values for R0
// in the graph construction, go ahead and define the
// realistic derived value as a pladeholder.
var R0 = (k * lambda * beta) / (a * d * u);
// We're using a range of values for R0. Generate a
// logarithmic range from 0.1 to 10 with 1.0 in the
// center. In addition to the R0 values, we also define
// a diverging color scale.
var R0Range = d3.range(-1,1.25,0.25).map(function(i) {
return Math.pow(10,i); }
);
var R0Colors = ['#006363','#2e8080','#619e9d','#98b9b9',
'#888','#d4a3b8','#ce749a','#bf4278','#a5004b'];
// Since we're going to graph the value of v on a
// logarithmic scale, define a minimum value greater
// than zero.
var vMin = 1e-8;
// Convenient ratios to convert from system state to
// meaningful physical quantities.
var ul2person = 5e6, // adults have about 5L of blood
t2days = 1 / d; // convert system time to days
// The system we're simulating is three differential
// equations. To keep the function definitions clean,
// we get the system parameter values from globals
// instead of function parameters.
function dx_dt(x,y,v) { return 1 - x - R0 * x * v; };
function dy_dt(x,y,v) { return R0 * x * v - alpha * y; };
function dv_dt(x,y,v) { return ( alpha * y - v ) / epsilon; };
// The time increment used in the simation.
var dt = 0.0001;
// Here's the function that actually simulates the system.
// It's inputs are the initial value for v, the number of
// points to generate, and the number of time steps between
// each point.
var simulate = function(v0, numPoints, numSteps) {
var result = [],
prev = { x: 1, y:0, v: v0, t: 0 };
for (var i=0; i<numPoints; i++) {
// At the beginning of each t-cycle, add a new
// data point to the results. Convert from the
// derived units used in the simulation to
// meaningful physical quantities when we do so.
result.push({
t: prev.t * t2days,
v: prev.v * ul2person
});
// Now quickly step through all the individual
// delta-t values for the current t-cycle.
for (var j=0; j<numSteps; j++) {
// Make sure none of the values go negative
// (and that v stays above the minimum).
var x = prev.x + dx_dt(prev.x, prev.y, prev.v) * dt,
y = prev.y + dy_dt(prev.x, prev.y, prev.v) * dt,
v = prev.v + dv_dt(prev.x, prev.y, prev.v) * dt;
// Make sure none of the values go negative
// (and that v stays above the minimum).
prev.x = Math.max(0,x);
prev.y = Math.max(0,y);
prev.v = Math.max(vMin,v);
prev.t += dt;
}
}
return result;
};
// Generate a separate dataset for each value of R0
// that we're graphing.
var datasets = R0Range.map(function(r0, i) {
// If we're less than half-way through the R0
// range, then start with a large virus population
// (a billion free virus particles). Otherwise,
// start with a single virus particle in the body.
var v0 = i <= (R0Range.length/2) ? 1e9/ul2person : 1/ul2person;
// Set the value for R0 and run the simulation.
R0 = r0;
return simulate(v0, width, 25);
});
// Define the scales that map a data value to a
// position on the vertical axis and to a time
// on the horizontal axis.
var v = d3.scale.log(),
t = d3.scale.linear();
// Set the scales
v.domain([d3.min(datasets, function(data) {
return d3.min(data, function(d) { return d.v; });
}), d3.max(datasets, function(data) {
return d3.max(data, function(d) { return d.v; });
})]).range([height,0]);
t.domain([0, datasets[0][datasets[0].length-1].t])
.range([0, width]);
// Define a convenience function to create a line on
// the chart. The horizontal axis (which D3 refers
// to as `x`) are the time values, and the vertical
// axis (which D3 refers to as `y`) are the x-values.
// The result is that `line` is function that, when
// passed a selection with an associated array of
// data points, returns an SVG path whose coordinates
// match the t- and x-scales of the chart.
var line = d3.svg.line()
.x(function(d) { return t(d.t); })
.y(function(d) { return v(d.v); });
// Create a selection for the data set
// and add the data to it. For each dataset,
// add a line to the chart.
var lines = g.selectAll(".line")
.data(datasets);
lines.enter().append("path")
.classed("line", true)
.attr("fill", "none")
.attr("stroke-width", "1px")
.attr("stroke", function(d,i) { return R0Colors[i]; })
.attr("d", line);
// Define a function that will create the axes
// when passed a data selection.
var vAxis = d3.svg.axis()
.scale(v)
.orient("left")
.ticks(5,"s");
var tAxis = d3.svg.axis()
.scale(t)
.tickSize(5, 0)
.orient("bottom");
// Add the v-axis.
svg.append("g")
.attr("class", "v axis")
.attr("transform", "translate("+margin.left+","+margin.top+")")
.call(vAxis)
.append("text")
.attr("x", -50)
.attr("y", -10)
.text("Free Virus Particles");
// Add the t-axis.
svg.append("g")
.attr("class", "t axis")
.attr("transform", "translate("+margin.left+","+(height+margin.top)+")")
.call(tAxis)
.append("text")
.attr("x", width/2)
.attr("y", 36)
.attr("text-anchor", "middle")
.text("Days Since Infection");
// Style the axis. We could do this with CSS, but
// by using JavaScript we ensure that the styles
// are embedded in the SVG itself. That's helpful
// if we ever want to extract the SVG from the web
// page as an image.
svg.selectAll(".axis line, .axis path")
.attr("fill", "none")
.attr("stroke", "#bbbbbb")
.attr("stroke-width", "2px")
.attr("shape-rendering", "crispEdges");
svg.selectAll(".axis text")
.attr("fill", "#444")
.attr("font-size", "14");
svg.selectAll(".axis .tick line")
.attr("stroke", "#d0d0d0")
.attr("stroke-width", "1");
// Finally, add the legends.
var legend = svg.append("g")
.attr("transform","matrix(1,0,0,1,60,428)");
legend.append("g")
.attr("transform","matrix(0.970129,0,0,1,-45.2062,-428)")
.append("rect")
.attr("x","533.172")
.attr("y","235.802")
.attr("width","95.9953")
.attr("height","184.219")
.attr("fill","none")
.attr("stroke-width","1.35")
.attr("stroke","#d0d0d0");
var label = legend.append("g")
.attr("transform","matrix(1,0,0,1,545.146,-174.78)");
label.append("text")
.attr("x","-67.361")
.attr("y","0")
.attr("font-size","14")
.attr("fill","#444444")
.text("R");
label.append("text")
.attr("x","-57.393")
.attr("y","4.662")
.attr("font-size","8.162")
.attr("fill","#444444")
.text("0");
label.append("text")
.attr("x","-52.2673")
.attr("y","0")
.attr("font-size","14")
.attr("fill","#444444")
.text(":");
legend.append("g")
.attr("transform","matrix(1,0,0,1,571.361,-174.78)")
.append("text")
.attr("x","-67.361")
.attr("y","0")
.attr("font-size","14")
.attr("fill","#444444")
.text("0.10");
legend.append("g")
.attr("transform","matrix(1,0,0,1,571.361,-74.78)")
.append("text")
.attr("x","-67.361")
.attr("y","0")
.attr("font-size","14")
.attr("fill","#444444")
.text("1.78");
legend.append("g")
.attr("transform","matrix(1,0,0,1,571.361,-154.78)")
.append("text")
.attr("x","-67.361")
.attr("y","0")
.attr("font-size","14")
.attr("fill","#444444")
.text("0.18");
legend.append("g")
.attr("transform","matrix(1,0,0,1,571.361,-54.78)")
.append("text")
.attr("x","-67.361")
.attr("y","0")
.attr("font-size","14")
.attr("fill","#444444")
.text("3.16");
legend.append("g")
.attr("transform","matrix(1,0,0,1,571.361,-134.78)")
.append("text")
.attr("x","-67.361")
.attr("y","0")
.attr("font-size","14")
.attr("fill","#444444")
.text("0.32");
legend.append("g")
.attr("transform","matrix(1,0,0,1,571.361,-34.78)")
.append("text")
.attr("x","-67.361")
.attr("y","0")
.attr("font-size","14")
.attr("fill","#444444")
.text("5.62");
legend.append("g")
.attr("transform","matrix(1,0,0,1,571.361,-114.78)")
.append("text")
.attr("x","-67.361")
.attr("y","0")
.attr("font-size","14")
.attr("fill","#444444")
.text("0.56");
legend.append("g")
.attr("transform","matrix(1,0,0,1,571.361,-94.78)")
.append("text")
.attr("x","-67.361")
.attr("y","0")
.attr("font-size","14")
.attr("fill","#444444")
.text("1.00");
legend.append("g")
.attr("transform","matrix(1,0,0,1,571.361,-14.78)")
.append("text")
.attr("x","-67.361")
.attr("y","0")
.attr("font-size","14")
.attr("fill","#444444")
.text("10.0");
legend.append("g")
.attr("transform","matrix(0.972202,-6.80854e-18,-6.80854e-18,1,19.1996,-460.396)")
.append("path")
.attr("d","M532.644,280.396l20.5719,0")
.attr("fill","none")
.attr("stroke-width","2.7")
.attr("stroke","#006363");
legend.append("g")
.attr("transform","matrix(0.972202,-6.80854e-18,-6.80854e-18,1,19.1629,-360.396)")
.append("path")
.attr("d","M532.644,280.396l20.5719,0")
.attr("fill","none")
.attr("stroke-width","2.7")
.attr("stroke","#d4a3b8");
legend.append("g")
.attr("transform","matrix(0.972202,-6.80854e-18,-6.80854e-18,1,19.1996,-440.396)")
.append("path")
.attr("d","M532.644,280.396l20.5719,0")
.attr("fill","none")
.attr("stroke-width","2.7")
.attr("stroke","#2e8080");
legend.append("g")
.attr("transform","matrix(0.972202,-6.80854e-18,-6.80854e-18,1,19.1629,-340.396)")
.append("path")
.attr("d","M532.644,280.396l20.5719,0")
.attr("fill","none")
.attr("stroke-width","2.7")
.attr("stroke","#ce749a");
legend.append("g")
.attr("transform","matrix(0.972202,-6.80854e-18,-6.80854e-18,1,19.1996,-420.396)")
.append("path")
.attr("d","M532.644,280.396l20.5719,0")
.attr("fill","none")
.attr("stroke-width","2.7")
.attr("stroke","#619e9d");
legend.append("g")
.attr("transform","matrix(0.972202,-6.80854e-18,-6.80854e-18,1,19.1629,-320.396)")
.append("path")
.attr("d","M532.644,280.396l20.5719,0")
.attr("fill","none")
.attr("stroke-width","2.7")
.attr("stroke","#bf4278");
legend.append("g")
.attr("transform","matrix(0.972202,-6.80854e-18,-6.80854e-18,1,19.1996,-400.396)")
.append("path")
.attr("d","M532.644,280.396l20.5719,0")
.attr("fill","none")
.attr("stroke-width","2.7")
.attr("stroke","#98b9b9");
legend.append("g")
.attr("transform","matrix(0.972202,-6.80854e-18,-6.80854e-18,1,19.1629,-300.396)")
.append("path")
.attr("d","M532.644,280.396l20.5719,0")
.attr("fill","none")
.attr("stroke-width","2.7")
.attr("stroke","#a5004b");
legend.append("g")
.attr("transform","matrix(0.972202,-6.80854e-18,-6.80854e-18,1,19.1996,-380.396)")
.append("path")
.attr("d","M532.644,280.396l20.5719,0")
.attr("fill","none")
.attr("stroke-width","2.7")
.attr("stroke","#888888");
</script>
</body>
</html>
Raw
thumbnail.lnk
http://jsdatav.is/img/thumbnails/hiv1.png
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment