# Innovative Approximate Analytical Solution for Standard Model of Viral Dynamics: Hepatitis C with Direct-Acting Agents as an Implemented Case

Authors: Hesham A. Elkaranshawy; Hossam M. Ezzat; Yasmine Abouelseoud; Nermeen N. Ibrahim
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

In this article, a novel approximate analytical solution is presented for solving the standard viral dynamic model. Basically, the standard model is used to study viral dynamics in patients for a wide range of viruses like HIV, HPV, HBV, and HCV. In this research work, the standard model for hepatitis C virus (HCV) is considered in detail; however, the analysis and results can be applicable for all other viruses. This standard model is used to study viral dynamics in patients treated with direct-acting antiviral agents (DAAs). Power series solution combined with Laplace–Padé resummation method (PSLP) is used to obtain the approximate analytical solutions for the model. To test the capability as well as the validity of the proposed method, results are compared with available viral load data published in the literature and with published simulation results. Good fits are obtained in the comparison for all cases considered. Given that medical specialists and physicians are more interested in solutions that yield direct and simple predictions, it is expected that the proposed approximate analytical solution would be attractive to them and help them to obtain a straightforward and a proper estimation regarding the viral load due to variations in treatment and/or patient’s parameters.


---

## Body

## 1. Introduction

Treatment for chronic hepatitis C infection began in the early 1990s with interferon-alpha [1]. This injectable drug worked by improving the immune system, rather than by specifically attacking the virus. In 1998, the oral drug ribavirin was added to interferon [2]. The development of the treatment occurred in 2002 with the approval of pegylated interferon-alpha, a process that makes interferon more durable and effective [3]. In 2011, antiviral medications that stop the virus reproducing, which called direct-acting antivirals (DAAs), appeared [4]. Sofosbuvir, commercial name Sovaldi, was approved in the United States in December 2013. Sovaldi-based oral therapy offers high cure rates for HCV infection with excellent tolerability [5]. Currently, the recommended DAA regimens are combinations of multidrug treatments, containing in particular NS5A, NS5B, NS4A inhibitors, such as glecaprevir + pibrentasvir, sofosbuvir and ribavirin (SOF + RBV), telaprevir + pegylated interferon and ribavirin (TVR + PR), sofosbuvir + simeprevir (SOF + SIM), and sofosbuvir + velpatasvir. For more details, see [6, 7, 8]. These regimens offer high cure rates with a rate of curation near 100%.The quick step of HCV drug development has led to the hopeful estimate that fully eradication of HCV is possible. Although there is still no vaccine for HCV, cure and eradication happened, but further study for the factors that increase this eradication rate is needed. The World Health Organization (WHO) has formulated the Global Health Sector Strategy on Viral Hepatitis, 2016–2021, and established service coverage targets to eliminate HCV as a public health threat by 2030 [9, 10]. Meanwhile, there remain many barriers that need to be overcome. Such barriers include the development of simplified and highly effective drug regimens, improving the rates of detection of infection, and the availability of funds including financial and medical expertise.Mathematical modeling is a useful tool in the study of virus dynamics for many types of viruses such as HCV, HBV, HPV, and HIV. It can be used to predict behavior under certain conditions or decide which parameters enhance the spread of disease. It may also be used to calculate the medications required to eradicate the disease or, at least, get it under control. These models can be used also to understand the biological mechanisms and interpret the experimental results [11]. Mathematical models in the form of a system of differential equations for the basic dynamics in vivo were developed and analyzed for HCV [8, 12–15], HBV [16, 17], HIV [18–22], and HPV [23]. Models for HCV treatment with DAAs therapy are considered in [8, 14, 15]. Authors in [24–26] developed a multiscale model that considered both intracellular viral RNA replication and extracellular viral infection.In general, numerical solutions were obtained for the models. However, analytical solutions could be useful for the estimation of parameters and for direct and simple predictions for the viral loads. Seldom, an analytical or an approximate analytical solution can be found in the literature. A simplified analytical solution for the standard model for HCV viral dynamics was constructed in [12]. It was assumed that the number of target cells is constant; hence, the system of ordinary differential equations (ODE) is reduced, linearized, and solved. It can be noticed that this system of ODE contains 7 parameters while the obtained solutions contain only 3 parameters. An analytical solution was obtained for the multiscale model in [24, 26]. All new infections during therapy were ignored in that analysis.In this research work, the approximate analytical solution is obtained for the standard model. Similar standard models have been used for HCV, HBV, and HIV in the literature [8, 12–23]. Nevertheless, the model considered is the HCV model. There is no reduction for the system of ODE or assumptions to simplify the equations of the model. Yet, the approximation is in the solution not in the system of differential equations. Therefore, it is more accurate and satisfactory than the previously mentioned analytical solutions. Power series solutions are first obtained for the system of ODE, and then, the Laplace–Padé resummation method (PSLP) is used to obtain the entailed approximate analytical solution. The technique used depends on the methodology developed in [27–34] for the general system of ODE. The solution is used for the analysis of the standard viral dynamic model of HCV after any type of DAA treatment initiation.


## 2. Standard Model

Consider the system of nonlinear ODEs for the standard viral dynamic mathematical model for HCV kinetics during treatment [11–14]:(1)dTdt=s−dT−βVT,(2)dIdt=βVT−δI,(3)dVdt=1−ερI−cV,where T is the target cells, which are produced at a constant rate s and are assumed to die per capita rate d, and are infected by virus, V, at rate β. Infected cells are assumed to die at per capita rate δ. Virions are generated at rate ρ per infected cell and cleared from serum at rate c per virion. Treatment is assumed to reduce the average viral production rate per infected cell from ρ to 1−ερ, where ε is the in vivo antiviral effectiveness of therapy (0<ε<1). At treatment initiation t=0, we used the standard assumption [10, 21, 22] that the system in the pretreatment is in steady state given by(4)dTdt=0,dIdt=0,dVdt=0,then,(5)T0=T0=cδβρ,I0=I0=−cdδ+βρsβδρ,V0=V0=−cdδ+βρsβcδ.


## 3. Solution Method

### 3.1. Power Series Solution

Assuming the following solutions of order N,(6)Tt=∑i=0NcTiti,(7)It=∑i=0NcIiti,(8)Vt=∑i=0NcViti.Substituting equations (6)–(8) into equations (1)–(3) and equating terms having the same powers of t. The cTi, cIi, and cVi coefficients can be calculated. For instance, the power series solution for Vt is(9)Vt=∑i=0NcViti=cV0+cV1t+cV2t2+cV3t3+⋯,where cV0=V0 is the initial condition of Vt and the coefficients cV1, cV2, and cV3 are as follows:(10)cV1=−cεc0,cV2=0.5c2εc0,cV3=16c0−c3ε−c2δε−cdδ+ρs−βc0cδ+ε2c2δ+cdδε−ρsβε+cδc0βε.


### 3.2. Laplace–Padé Resummation

The Padé approximant for a series function is a function with the same power series expansions as the original series function. It is suitable for approximating a divergent series function. The approximant is derived by expanding the function as a ratio of two power series and determining both the numerator and denominator coefficients [28–34]. If we have a function fu that can be represented in a power series form as(11)fu=∑i=0∞cuiui.A Padé approximant is a rational function:(12)LM=b0+b1u+b2u2+⋯+bLuLd0+d1u+d2u2+⋯+dMuM,which has a Maclaurin expansion that agrees with the series of equation (11) approximately.d0 is chosen to be equal to 1. So there are L+1 independent numerator coefficients and M independent denominator coefficients, making L+M+1 unknown coefficients. The L/M should fit the power series of equation (11) through the orders 1, u, u2, … , uL+M. Hence,(13)∑i=0∞cuiui=b0+b1u+b2u2+⋯+bLuL1+d1u+d2u2+⋯+dMuM+OuL+M+1.Up to the term of orderL+M, we can write(14)cu0+cu1u+⋯+cuL+MuL+M1+d1u+d2u2+⋯+dMuM=b0+b1u+b2u2+⋯+bLuL.Denominator coefficientsd1,d1,…,dM can be found from equation (14) by equating the coefficients of uL,uL+1,…,uL+M. Numerator coefficients b0+b1+⋯bL can be found from equation (14) by equating the coefficients of u0,u,u2,…,uL. Hence, we have constructed the L/M Padé approximant which agrees with ∑i=0∞cuiui through the order uL+M.To extend the convergence zone of the power series solution, Laplace–Padé resummation can be introduced. In Laplace–Padé resummation, the Laplace transform is applied to the given series function, then Padé approximate is obtained for the resulted function, and consequently, the inverse Laplace transform is used to specify the final solution in the form of exponential functions. To apply Laplace–Padé resummation to the series solution of the viral loadVt, the following procedure is employed:(1)
From equation (9), the third-order power series solution of Vt is(15)Vt=cV0+cV1t+cV2t2+cV3t3.(2)
Taking Laplace transform to the series solution in equation (15) to obtain(16)Fs=cV0s+cV1s2+2cV2s3+6cV3s4.(3)
1/u is written instead of s in equation (16) to obtain(17)Zu=cV0u+cV1u2+2cV2u3+6cV3u4.(4)
2/2 Padé approximant for equation (17) can be written as(18)2/2=b0+b1u+b2u21+d1u+d2u2.Then,(19)cV0u+cV1u2+2cV2u3+6cV3u4=b0+b1u+b2u21+d1u+d2u2.Then, equation (19) can be written as(20)cV0u+cV1u2+2cV2u3+6cV3u41+d1u+d2u2=b0+b1u+b2u2.Denominator coefficientsd1 and d2 can be found by equating the coefficients of u3 and u4. Numerator coefficients b0,b1, and b2 can be found by equating the coefficients of u0,u, and u2. Using Mathematica program, hence(21)b0=0,b1=cV0,b2=cV13−4cV1cV2cV0+6cV3cV02f,d1=−2cV1cV2−3cV3cV0cV12−2cV2cV0,d2=−2−2cV22+3cV1cV3cV12−2cV2cV0.(5)
1/s is written instead of u in 2/2; hence, Fs is obtained as(22)Fs=cV0/s+b2/s21+d1/s+d2/s2.(6)
Using inverse Laplace transform forFs, the viral load Vt is obtained as(23)Vt=A1e−D1∗t+cV0−A1e−D2∗t,where(24)A1=12B−cV13+3cV1cV2cV0−3cV3cV02+cV0B,D1=−A−Bf,D2=−A+Bf,f=cV12−2cV2cV0,A=cV1cV2−3cV3cV0,B=−3cV12cV22+6cV13cV3+8cV23cV0−18cV1cV2cV3cV0+9cV32cV02.Equation (23) combined with equations (10) and (24) represents the general approximate analytical solution for the viral load for the standard dynamic model given in equations (1)–(3). It is worth to notice that this solution contains the 7 parameters given in the system of equations (1)–(3). Hence, it reflects properly the biological features inherited in the system. It can be noticed that the simplified approximate solution given in equation (23)  contains only 3 parameters.


## 3.1. Power Series Solution

Assuming the following solutions of order N,(6)Tt=∑i=0NcTiti,(7)It=∑i=0NcIiti,(8)Vt=∑i=0NcViti.Substituting equations (6)–(8) into equations (1)–(3) and equating terms having the same powers of t. The cTi, cIi, and cVi coefficients can be calculated. For instance, the power series solution for Vt is(9)Vt=∑i=0NcViti=cV0+cV1t+cV2t2+cV3t3+⋯,where cV0=V0 is the initial condition of Vt and the coefficients cV1, cV2, and cV3 are as follows:(10)cV1=−cεc0,cV2=0.5c2εc0,cV3=16c0−c3ε−c2δε−cdδ+ρs−βc0cδ+ε2c2δ+cdδε−ρsβε+cδc0βε.


## 3.2. Laplace–Padé Resummation

The Padé approximant for a series function is a function with the same power series expansions as the original series function. It is suitable for approximating a divergent series function. The approximant is derived by expanding the function as a ratio of two power series and determining both the numerator and denominator coefficients [28–34]. If we have a function fu that can be represented in a power series form as(11)fu=∑i=0∞cuiui.A Padé approximant is a rational function:(12)LM=b0+b1u+b2u2+⋯+bLuLd0+d1u+d2u2+⋯+dMuM,which has a Maclaurin expansion that agrees with the series of equation (11) approximately.d0 is chosen to be equal to 1. So there are L+1 independent numerator coefficients and M independent denominator coefficients, making L+M+1 unknown coefficients. The L/M should fit the power series of equation (11) through the orders 1, u, u2, … , uL+M. Hence,(13)∑i=0∞cuiui=b0+b1u+b2u2+⋯+bLuL1+d1u+d2u2+⋯+dMuM+OuL+M+1.Up to the term of orderL+M, we can write(14)cu0+cu1u+⋯+cuL+MuL+M1+d1u+d2u2+⋯+dMuM=b0+b1u+b2u2+⋯+bLuL.Denominator coefficientsd1,d1,…,dM can be found from equation (14) by equating the coefficients of uL,uL+1,…,uL+M. Numerator coefficients b0+b1+⋯bL can be found from equation (14) by equating the coefficients of u0,u,u2,…,uL. Hence, we have constructed the L/M Padé approximant which agrees with ∑i=0∞cuiui through the order uL+M.To extend the convergence zone of the power series solution, Laplace–Padé resummation can be introduced. In Laplace–Padé resummation, the Laplace transform is applied to the given series function, then Padé approximate is obtained for the resulted function, and consequently, the inverse Laplace transform is used to specify the final solution in the form of exponential functions. To apply Laplace–Padé resummation to the series solution of the viral loadVt, the following procedure is employed:(1)
From equation (9), the third-order power series solution of Vt is(15)Vt=cV0+cV1t+cV2t2+cV3t3.(2)
Taking Laplace transform to the series solution in equation (15) to obtain(16)Fs=cV0s+cV1s2+2cV2s3+6cV3s4.(3)
1/u is written instead of s in equation (16) to obtain(17)Zu=cV0u+cV1u2+2cV2u3+6cV3u4.(4)
2/2 Padé approximant for equation (17) can be written as(18)2/2=b0+b1u+b2u21+d1u+d2u2.Then,(19)cV0u+cV1u2+2cV2u3+6cV3u4=b0+b1u+b2u21+d1u+d2u2.Then, equation (19) can be written as(20)cV0u+cV1u2+2cV2u3+6cV3u41+d1u+d2u2=b0+b1u+b2u2.Denominator coefficientsd1 and d2 can be found by equating the coefficients of u3 and u4. Numerator coefficients b0,b1, and b2 can be found by equating the coefficients of u0,u, and u2. Using Mathematica program, hence(21)b0=0,b1=cV0,b2=cV13−4cV1cV2cV0+6cV3cV02f,d1=−2cV1cV2−3cV3cV0cV12−2cV2cV0,d2=−2−2cV22+3cV1cV3cV12−2cV2cV0.(5)
1/s is written instead of u in 2/2; hence, Fs is obtained as(22)Fs=cV0/s+b2/s21+d1/s+d2/s2.(6)
Using inverse Laplace transform forFs, the viral load Vt is obtained as(23)Vt=A1e−D1∗t+cV0−A1e−D2∗t,where(24)A1=12B−cV13+3cV1cV2cV0−3cV3cV02+cV0B,D1=−A−Bf,D2=−A+Bf,f=cV12−2cV2cV0,A=cV1cV2−3cV3cV0,B=−3cV12cV22+6cV13cV3+8cV23cV0−18cV1cV2cV3cV0+9cV32cV02.Equation (23) combined with equations (10) and (24) represents the general approximate analytical solution for the viral load for the standard dynamic model given in equations (1)–(3). It is worth to notice that this solution contains the 7 parameters given in the system of equations (1)–(3). Hence, it reflects properly the biological features inherited in the system. It can be noticed that the simplified approximate solution given in equation (23)  contains only 3 parameters.


## 4. Study Cases

To illustrate the capabilities and competences of the proposed method, four cases of study are presented. These cases consider the analysis of viral kinetics using the standard viral dynamic mathematical model of HCV. The proposed PSLP method is applied to solve the nonlinear dynamic model of the viral kinetics for some patients after the initiation of treatment with DAAs.For the first two study cases, the predictions of the PSLP method are compared, for each patient, with viral load data available in the literature [24, 25]. To provide best fits of data for each patient, parameters δ,ε,c, and ρ are estimated using the relation for initial condition V0=−cdδ+βρs/βcδ and equation (23). In the first case, patients were infected with HCV genotype 1 and treated with danoprevir. The viral load for these patients is checked for 13 days after initiation of danoprevir, the data are available in [24], and the parameters are given in Table 1. In the second case, patients were treated with daclatasvir and the viral load of patients is checked for 2 days after initiation of daclatasvir, the data are available in [25], and the parameters are given in Table 2. Figures 1 and 2 demonstrate the comparison between the solution of the PSLP method and the corresponding viral load data for each patient.Table 1
Parameter values used for the patients treated with danoprevir.
Patient
c
V
0
log
10
IU
/
mL
c
day
−
1
ε
δ
day
−
1
ρ
day
−
1
01–94GK
7.42
7.38
0.9995
0.15
149.684
03–94HD
6.72
12.44
0.998
0.29
151.188
03–94EA
5.79
10.5
0.998
0.17
11.212
03–94KG
6.98
9.4
0.98
0.35
246.748
04–94XD
6.63
10.26
0.9995
0.33
116.31Table 2
Parameter values used for the patients treated with daclatasvir.
Patient
c
V
0
log
10
IU
/
mL
c
day
−
1
ε
δ
day
−
1
ρ
day
−
1
PAT 8
5.64
31.3
0.996
1.89
289.649
PAT 42
5.65
21
0.998
0.39
40.741
PAT 68
7.15
20.2
0.998
1.27
2826.95
PAT 69
6.14
22
0.999
0.75
278.629
PAT 83
5.45
21.8
0.996
1.01
81.608Figure 1
Comparison between the approximate analytical solution using the PSLP method and viral data for patients treated with danoprevir.s=1.3∗105 cells/ml, d=0.01day−1, and β=5∗10−8mlday−1virion−1 and the rest of the parameter values are given in Table 1. (a) 01–94GK, (b) 03–94HD, (c) 03–94EA, (d) 03–94KG, and (e) 04–94XD.
(a)
(b)
(c)
(d)
(e)Figure 2
Comparison between the approximate analytical solution using the PSLP method and viral data for patients treated with daclatasvir.s=1.3∗105 cells/ml, d=0.01day−1, and β=5∗10−8mlday−1virion−1 and the rest of the parameter values are given in Table 2. (a) PAT 8, (b) PAT 42, (c) PAT 68, (d) PAT 69, and (e) PAT 83.
(a)
(b)
(c)
(d)
(e)Two other cases representing recently used combination of multidrug DAA treatments are considered. The predictions of the PSLP method are compared with the corresponding viral load obtained by simulation in [7]. The used parameter values are given in Table 3, where δ,εandc have the same values assigned in [7] and ρ has been chosen accordingly. The PSLP results and the published simulation results for 25 days after the initiation of treatment with TVR + PR are shown in Figure 3. The PSLP results and the published simulation results for 28 days after the initiation of treatment with SOF + SIM are shown in Figure 4. In both cases, the PSLP solution and the simulation results are almost identical. The comparison shows that the proposed PSLP method provides adequate approximate analytical solutions using the standard viral dynamic model for all the considered cases. It is worth to mention that, in these two cases, the rate of curation is near 100%.Table 3
Parameter values used with SOF + SIM and TVR + PR.
DAAs
c
V
0
log
10
IU
/
mL
c
day
−
1
ε
δ
day
−
1
ρ
day
−
1
SOF + SIM
5.98
5.28
0.998
0.27
8.18
TVR + PR
5.98
5.28
0.999
0.27
8.18Figure 3
Comparison between the results obtained by the PSLP method and by simulation in [7] for treatment with TVR + PR. s=1.3∗105 cells/ml, d=0.01day−1, and β=5∗10−8mlday−1virion−1 and the rest of the parameter values are given in Table 3.Figure 4
Comparison between the results obtained by the PSLP method and by simulation in [7] for treatment with SOF + SIM. s=1.3∗105 cells/ml, d=0.01day−1, and β=5∗10−8mlday−1virion−1 and the rest of the parameter values are given in Table 3.Since it is highly desirable to predict the patient’s response to a specific treatment regimen prior to the starting of the treatment itself, the PSLP solution offers a simple and powerful tool for medical specialists and physicians to perform this task. They can use the patient’s parameters to calculate the constants in equations (10) and (24) and substitute them in equation (23) to get a closed-form solution for the viral load. Hence, the viral load can be plotted versus time or simply the viral load can be estimated at any instant by direct substitution in the closed-form solution.


## 5. Conclusions

Power series solution combined with the Laplace–Padé resummation method (PSLP) has been used to obtain a general approximate analytical solution for the nonlinear standard viral dynamic model of HCV for patients treated with DAAs. However, the solution is not limited to HCV model, and it can be applied to other viruses like HIV and HBV, for example, given that proper parameters are used. To test the applicability and accuracy of the proposed method, results have been compared with viral load data and with published simulated results. Satisfactory agreement between the PSLP solution and the corresponding viral load data has been found for all the considered cases. The PSLP results and the published simulated results are almost identical. The comparison proves that this innovative PSLP solution can be used with confidence for solving the nonlinear standard viral dynamic model. This solution can conveniently be used to fit patient data and estimate parameter values. So it would facilitate for physicians to monitor the changes in the viral load due to changes in treatment and to deal with changes in patient’s parameters.



---

*Source: 1454739-2019-09-16.xml*