September 22, 2017, Friday

# MBW:Fundamental Mathematics of the Delay Differential Equations in the HIV-1 Virus Model

This is an article is written by Sai'd Azzam and John Nelson Kane

## Executive Summary

In this mathematical model we attempt to better understand the dynamics of the HIV-1 infection by incorporating an intracellular delay, which is the time period between an antiviral drug reaching a cell and the effectiveness of the drug. In the past, models have made the unrealistic assumption that the intracellular delay is constant or there is no delay after ingestion of the drug. This model introduces this delay in the form of a continuous distribution with the underlying principle that different cells will be associated with a different delay time. One fundamental equation of this model is known as an integro-differential equation. The focus of this study will be on the analytic methods used to transform this integro-differential equation into a set of linear 1st order differential equations in which analytical solutions can be obtained. This will be done with a method known as the linear chain trick. By incorporating this method it will be shown that a time delay equation can be transformed into an ordinary equation of only current time, thus making analytic solutions possible. As a final conclusion we will present the final set of 1st order linear differential equations that will equally describe this model.

## Introduction

The mathematics describing this model are motivated by the research and attempted understanding of the Human Immunodeficiency Virus (HIV). HIV is a replicating virus and the cause of acquired immunodeficiency syndrome (AIDS). HIV-1 is the most common and pathogenic strain of the virus and is typically associated with three stages. Once a patient is first infected, the amount of virus found in the blood stream is relatively high. This high concentration may last several weeks to a few months at which point the concentration decreases to some stable lower level. Following the beginning of this stable level of virus concentration, a typically long period follows where the CD4+ T-cell count slowly decreases. A CD4+ T-cell is also known as a T helper cell that is a type of white blood cell used to help the activity of other immune cells. This usually long period of a decrease in T-cells appears to be a quasi-steady state, because although it appears to be stable, the virus concentration is noticed to change slightly when analyzed on a long time scale. This quasi-steady state can last for up to a few years with little change. In the final stage of the disease, the CD4+ T-cell count drops to a point in which the immune system fails and the viral load increases rapidly, usually causing fatality. Much research has been done to better understand this middle period in which the T-cells slowly decrease and the virus concentration stays at a relatively stable value. If one can provide the correct combination of antiretroviral drugs with treatment to better stabilize the white blood cells than one could dramatically prolong life and even essentially make the disease nonfatal to the human body. This model attempts to analytically describe the middle stage, the quasi-steady state of the virus. Previous models have attempted to describe this stage with either a constant delay, or the special case of no delay at all. These are not realistic for obvious reasons, although constant delays allow for analytic solutions to be easily obtain for the differential equations in those models are first order and linear. This model will incorporate a gamma distribution for the time delay with the physical interpretation that there is no constant time for which the decay occurs in any given cell. The density function of a gamma distribution is shown below:

To account for all possible times of the delay one must integrate over the entire gamma distribution, thus an integrand will appear in the new model. This part of the new model is known as an integro-differential equation. The abstract mathematical methods used to interpret this integro-differential equation will be the core of this paper. We will apply a method known as the linear chain trick to this integro-differential equation to transform it into equations that are more easily understood and more convenient to work with numerically.

## Delay / Background

There are two different types of delays that will effect the infectious/treatment process: a pharmacological delay and an intracellular delay. A pharmacological delay is the time between the ingestion of the drug and it’s first appearance within the cells. An intracellular delay is the time between the initial infection of a cell and the release of new viruses. Studies show from Perelson et al. that the pharmacological delay is around two hours while the intracellular delay is of the magnitude of a few days. It is also must be considered that the pharmacological delay only occurs at the commencement of treatment, while intracellular delay occurs continuously. For this reason we can ignore pharmacological delay and we will do so for the following model.

The intracellular delay will be in the form of a gamma distribution. As shown in the following methods, this will allow for an analytical solution to be found. The gamma distribution is also considered to be a reasonable delay distribution within a cell from the biological sense.

## Mathematical Model With Variable/Parameter Definitions

The main goal of this paper is to have a biologically plausible answer, and to analytically transform an integro-differential equation into a system of first order differential equations. The parameters associated with the model are: the density of productively infected cells, I, the concentration of infections virus, V(I), the viral clearance rate constant per day, c, and the concentration of virus has been rendered non-infections by the protease inhibitors, V(NI).

The total concentration of virus is V=V(I) + V(NI). We can also assume that the density of the non-infected cells, T, will remain constant. This is a fair assumption since the data are only examined for a short period of time following the initial ingestion of the drug.

In the absence of delay, Perelson et al. proposed the following model:

$(1)\qquad \qquad {\frac {dI}{dt}}=kTV_{I}-\delta I$

$(2)\qquad \qquad {\frac {dV_{I}}{dt}}=(1-\eta )pI(t)-cV_{I}(t)$

$(3)\qquad \qquad {\frac {V_{{NI}}}{dt}}=\eta pI(t)-cV_{{NI}},$

We make note that k is the infection rate constant, p is the rate at which a productively infected cell release viruses. Productivity infected cells die at a rate $\delta$ per cell. Plasma viruses are cleared at rate c per virus.$\eta$ respresents the drug efficiency. If $\eta =0$ than the drug has no effective and if $\eta =1$ than the drug is 100% effective.

As noted above this model will only account for the intracellular delay, and the pharmacological delay will be ignored. With the addition of the delay the following model is proposed in the form of integro-differential equations as follows:

$f(t')=g_{{n,b}}(t')={\frac {t'^{{(n-1)}}}{(n-1)!b^{n}}}e^{{\frac {-t'}{b}}}$

We multiplied the gamma distribution, f(t'), by the change of the infectious virus equation, ${\frac {dI}{dt}}$, and integrate to find the probability.

$(4a)\qquad \qquad {\frac {dI}{dt}}=\int _{0}^{\infty }kTf(t')V_{I}(t-t')e^{{-mt'}}dt-\delta I(t),$

The equations (4b) and (4c) are already in the form of first order differential equations, thus we will not transform them.

$(4b)\qquad \qquad {\frac {dV_{I}}{dt}}=[1-h(t)](1-\eta )pI(t)-cV_{I}(t),$

$(4c)\qquad \qquad {\frac {V_{{NI}}}{dt}}=h(t)\eta pI(t)-cV_{{NI}},$

Throughout the analysis we will assume that the system being treated is at steady state prior to being treated. We will also make note of $t=0$ at the time the drugs initially takes effect. This allows for the ignorance of the pharmacological delay. Noting that there are zero non-infectious virus cells at $t=0$ we can allow for initial conditions as follows:

$(5)\qquad \qquad I(0)={\frac {c}{p}}V_{I}(0)$

$(6)\qquad \qquad T=T(0)={\frac {c\delta }{kp}}$

The intracellular delay is not always known to an exact, nor is it the same value for every cell undergoing the therapy process. In fact it is unrealistic to represent this delay with a constant thus is will be represented with a gamma distribution known as the delay distribution .

with a mean, $\tau =nb$, a variance, $nb^{2}$ and a max peak, $(n-1)b$. This distribution is widely used in mathematical biology because it is a plausible delay and allows for analytic solutions to be obtained.

### Converting time delay defferential equations into ordinary differential equation

From 4(a) we wish to eliminate the term $e^{{-mt'}}$ , as a first step to approach an analytical solution, thus we make a substitution. We will define two new variables ${\hat {b}}$ and ${\hat {k}}$

$(8)\qquad \qquad {\hat {b}}\equiv {\frac {b}{1+mb}}$ and

${\hat {k}}\equiv {\frac k{(1+mb)^{n}}}$

We need to solve for b in terms of ${\hat {b}}$ and ${\hat {k}}$ thus,

$b={\hat {b}}(1+mb)={\hat {b}}+mb{\hat {b}}$

$b(1-m{\hat {b}})={\hat {b}}$

$b={\frac {{\hat {b}}}{(1-m{\hat {b}})}}$

Here we will solve for k , in term of ${\hat {b}}$ and ${\hat {k}}$

$k={\hat {k}}(1+mb)^{n}$

$(1+mb)={\frac {b}{{\hat {b}}}}={\frac {1}{1-m{\hat {b}}}}$

$k={\hat {k}}{\frac {1}{1-m{\hat {b}}}}$

Our integro-differential equation is of the form:

${\frac {dI}{dt}}=kT\int _{0}^{\infty }{\frac {t'^{{(n-1)}}}{(n-1)!b^{n}}}e^{{\frac {-t'}{b}}}V_{I}(t-t')e^{{-mt'}}dt'-\delta I(t)$

Substituting for b and k with the expressions above in terms of ${\hat {b}}$ and ${\hat {k}}$, and we get:

${\frac {dI}{dt}}={\hat {k}}T\int _{0}^{\infty }{\frac {(1-m{\hat {b}})^{n}t'^{{(n-1)}}}{(n-1)!{\hat {b}}^{n}(1-m{\hat {b}})^{n}}}e^{{\frac {-t'}{{\hat {b}}}}}e^{{mt'}}V_{I}(t-t')e^{{-mt'}}dt'-\delta I(t)$

As we can see, we cancel the $e^{{-mt}}$ and $(1-m{\hat {b}})^{n}$ terms. Thus we arrive at:

${\frac {dI}{dt}}={\hat {k}}T\int _{0}^{\infty }{\frac {t'^{{(n-1)}}}{(n-1)!({\hat {b}})^{n}}}e^{{\frac {-t'}{{\hat {b}}}}}V_{I}(t-t')dt'-\delta I(t)$

Thus our equation is now in terms of ${\hat {b}}$ and ${\hat {k}}$. The equation is now similar to (4a) with the exception of ${\hat {b}}$ in replacement of b, ${\hat {k}}$ in replacement of k and no $e^{{-mt}}$ term. We now have a gamma distribution in terms of ${\hat {b}}$ and is defined as:

$g_{{n,{\hat {b}}}}(t')={\frac {t'^{{(n-1)}}}{(n-1)!{\hat {b}}^{n}}}e^{{\frac {-t'}{{\hat {b}}}}}$

We now introduce the function $E_{n}(t)$.

This $E_{n}(t)$ term is defined as:

$E_{n}(t)=\int _{0}^{\infty }g_{{n,{\hat {b}}}}(t')V_{I}(t-t')dt'$

With the gamma distribution expanded we have:

$E_{n}(t)=\int _{0}^{\infty }{\frac {t'^{{(n-1)}}}{(n-1)!{\hat {b}}^{n}}}e^{{\frac {-t'}{{\hat {b}}}}}V_{I}(t-t')dt'$

To begin the linear chain trick, we introduce a new variable u = t-t' and substitute it into $E_{n}(t)$

$u=t-t'$

solving for t', we get:

$t'=u-t$

After substituting u into $E_{n}(t)$ we get:

$E_{n}(t)=\int _{{-\infty }}^{t}{\frac {(t-u)^{{n-1}}}{(n-1)!{\hat {b}}^{n}}}e^{{{\frac {-(t-u)}{{\hat {b}}}}}}V_{I}(u)du$

Then we differentiate the integrand with respect to the parameter t. We notice that the integral is over du thus we can use Leibniz's rule of differentiating integrals. This derivative with respect to t, applying the chain rule is:

${\dot {E}}_{n}(t)=\int _{{-\infty }}^{t}{\frac {(n-1)(t-u)^{{n-2}}}{(n-1)!{\hat {b}}^{n}}}e^{{{\frac {-(t-u)}{{\hat {b}}}}}}V_{I}(u)du$

$\qquad -\int _{{-\infty }}^{t}{\frac {(t-u)^{{n-1}}}{(n-1)!{\hat {b}}^{n}}}{\frac {e^{{{\frac {-(t-u)}{{\hat {b}}}}}}}{{\hat {b}}}}V_{I}(u)du$

$\qquad +{\frac {(t-t)^{{n-1}}}{(n-1)!{\hat {b}}^{n}}}e^{{{\frac {-(t-t)}{{\hat {b}}}}}}V_{I}(t)$

Then we simplify the terms and we get:

${\dot {E}}_{n}(t)=\int _{{-\infty }}^{t}{\frac {(t-u)^{{n-2}}}{(n-2)!{\hat {b}}^{n}}}e^{{{\frac {-(t-u)}{{\hat {b}}}}}}V_{I}(u)du$

$\qquad -\int _{{-\infty }}^{t}{\frac {(t-u)^{{n-1}}}{(n-1)!{\hat {b}}^{{n+1}}}}e^{{{\frac {-(t-u)}{{\hat {b}}}}}}V_{I}(u)du+g_{{n,{\hat {b}}}}(0)V_{I}(t)$

Since $g_{{n,{\hat {b}}}}(0)=0$ , we get

${\dot {E}}_{n}(t)=\int _{{-\infty }}^{t}{\frac {(t-u)^{{n-2}}}{(n-2)!{\hat {b}}^{n}}}e^{{{\frac {-(t-u)}{{\hat {b}}}}}}V_{I}(u)du$

$\qquad -\int _{{-\infty }}^{t}{\frac {(t-u)^{{n-1}}}{(n-1)!{\hat {b}}^{{n+1}}}}e^{{{\frac {-(t-u)}{{\hat {b}}}}}}V_{I}(u)du$

If we look closely, we notice that the above terms mean,

${\dot {E}}_{n}(t)={\frac {\left[E_{{n-1}}(t)-E_{n}(t)\right]}{{\hat {b}}}}$

This is only valid for $n>1$ thus we have to get a solution in term of $n=1$. Thus we plug 1 for the $E_{n}(t)$ equation, and we get:

$E_{1}(t)=\int _{{-\infty }}^{t}{\frac {1}{{\hat {b}}}}e^{{{\frac {-(t-u)}{{\hat {b}}}}}}V_{I}(u)du$

We differentiate in the same manner with respect to the parameter t. First we manipulate for simplification. We pull the term $e^{{{\frac {-t}{{\hat {b}}}}}}$ out of the intergal, which is considered a constant in terms of the variable u.

$E_{1}(t)=e^{{{\frac {-t}{{\hat {b}}}}}}\int _{{-\infty }}^{t}{\frac {1}{{\hat {b}}}}e^{{{\frac {u}{{\hat {b}}}}}}V_{I}(u)du$

Using Leibniz's rule of differentiation of an integrand we differentiate with respect to t.

${\dot {E}}_{1}(t)={\frac {e^{{{\frac {-t}{{\hat {b}}}}}}}{{\hat {b}}}}\int _{{-\infty }}^{t}{\frac {1}{{\hat {b}}}}e^{{{\frac {u}{{\hat {b}}}}}}V_{I}(u)du+e^{{{\frac {-t}{{\hat {b}}}}}}{\frac {1}{{\hat {b}}}}e^{{{\frac {t}{{\hat {b}}}}}}V_{I}(t)$

${\dot {E}}_{1}(t)={\frac {V_{I}(t)}{{\hat {b}}}}-{\frac {1}{{\hat {b}}}}\int _{{-\infty }}^{t}e^{{\frac {-(t-u)}{{\hat {b}}}}}V_{I}(u)du$

Again, if we look closely ${\frac {dE_{1}}{dt}}$ is just:

${\frac {dE_{1}}{dt}}={\frac {\left[V_{I}(t)-E_{1}(T)\right]}{{\hat {b}}}}$

$(9)\qquad \qquad {\frac {dI}{dt}}={\hat {k}}T\int _{0}^{\infty }g_{{n,{\hat {b}}}}(t')V_{I}(t-t')dt'-\delta I(t)$

After all of the above manipulations, now the integro-differential equations given in (4a), (4b) and (4c) are now equivalent to the following system of ordinary 1st order linear differential equations.

$(10a)\qquad \qquad {\frac {dI}{dt}}={\hat {k}}TE_{n}(t)-\delta I(t)$

$(10b)\qquad \qquad {\frac {dE_{1}}{dt}}={\frac {[V_{I}(t)-E_{1}(t)]}{{\hat {b}}}},$

$(10c)\qquad \qquad {\frac {dE_{i}}{dt}}={\frac {[E_{{i-1}}(t)-E_{i}(t)]}{{\hat {b}}}},\qquad i=2,......,n$

From (4b), our initial linear differential equation is:

${\frac {dV_{I}}{dt}}=[1-h(t)](1-\eta )pI(t)-cV_{I}(t)$

We will assume 100% effectiveness ( $\eta =1$), thus we arrive at:

$(10d)\qquad \qquad {\frac {dV_{I}}{dt}}=-cV_{I}(t)$

From (4c), our initial linear differential equation is:

${\frac {V_{{NI}}}{dt}}=h(t)\eta pI(t)-cV_{{NI}}(t),$

Assuming 100% effectiveness ( $\eta =1$), we get:

$(10e)\qquad \qquad {\frac {V_{{NI}}}{dt}}=h(t)pI(t)-cV_{{NI}}(t),$

$Ei(t)\qquad$ for $i=1,......,n$ are auxiliary functions and are related to the delay distribution as follows:

$(11)\qquad \qquad E_{i}(t)={\hat {k}}T\int _{0}^{\infty }g_{{i,{\hat {b}}}}(t')V_{I}(t-t')dt'$

## Results and Discussion

We began with a system of three delay differential equations, one being an integro-differential equation. Throughout this paper we have shown that through a series of steps and a method known as the linear chain trick, the integro-differential equation can be re-written in terms of a system of 1st order differential equations that can be solved for analytically. Our delay is in the form of a gamma distribution. One must account for all possibilities of the delay, thus it is needed to integrate over the whole distribution density function. This is the reason for an integrand appearing in the model to give us an integro-differential equation. The initial model is shown below with the integro-differential equation being the first presented.

$\qquad \qquad {\frac {dI}{dt}}=\int _{0}^{\infty }kTf(t')V_{I}(t-t')e^{{-mt'}}dt-\delta I(t),$

$\qquad \qquad {\frac {dV_{I}}{dt}}=[1-h(t)](1-\eta )pI(t)-cV_{I}(t),$

$\qquad \qquad {\frac {V_{{NI}}}{dt}}=h(t)\eta pI(t)-cV_{{NI}},$

It has been shown through that with the linear chain trick, the integro-differential equation can be transformed into a system of 1st order differential equations. Noticing that the original 2nd and 3rd equations in model are already in linear form they will remain in the final model. The only difference is that the drug has been assumed to be 100% effective. It has been shown that the above model presented is equivalent to the system of 1st order linear differential equations shown below:

$\qquad \qquad {\frac {dI}{dt}}={\hat {k}}TE_{n}(t)-\delta I(t)$

$\qquad \qquad {\frac {dE_{1}}{dt}}={\frac {[V_{I}(t)-E_{1}(t)]}{{\hat {b}}}},$

$\qquad \qquad {\frac {dE_{i}}{dt}}={\frac {[E_{{i-1}}(t)-E_{i}(t)]}{{\hat {b}}}},\qquad i=2,......,n$

$\qquad \qquad {\frac {dV_{I}}{dt}}=-cV_{I}(t)$

$\qquad \qquad {\frac {V_{{NI}}}{dt}}=h(t)pI(t)-cV_{{NI}}(t),$

Now that a first order system has been obtained, one can go further by working towards an analytic result. Noticing that there are now n + 3 linear equations one might find it more useful to solve using a software package with the n that is desired. To further investigate how the mathematics presented here will model a biological system such as HIV-1 dynamics and the associated delay, one could study the effects of given parameters with a sensitivity analysis.