Homoclinic and heteroclinic solutions to a hepatitis C evolution model

Tadas Telksnys; Zenonas Navickas; Romas Marcinkevicius; Maosen Cao; Minvydas Ragulskis

doi:10.1515/math-2018-0130

Enjoy 40% off

academic books on De Gruyter Brill *

Article Open Access

Homoclinic and heteroclinic solutions to a hepatitis C evolution model

Tadas Telksnys , Zenonas Navickas , Romas Marcinkevicius , Maosen Cao and Minvydas Ragulskis

Published/Copyright: December 31, 2018

Published by

Become an author with De Gruyter Brill

Submit Manuscript Author Information Explore this Subject

$Open Mathematics$

From the journal Open Mathematics Volume 16 Issue 1

Abstract

Homoclinic and heteroclinic solutions to a standard hepatitis C virus (HCV) evolution model described by T. C. Reluga, H. Dahari and A. S. Perelson, (SIAM J. Appl. Math., 69 (2009), pp. 999–1023) are considered in this paper. Inverse balancing and generalized differential techniques enable derivation of necessary and sufficient existence conditions for homoclinic/heteroclinic solutions in the considered system. It is shown that homoclinic/heteroclinic solutions do appear when the considered system describes biologically significant evolution. Furthermore, it is demonstrated that the hepatitis C virus evolution model is structurally stable in the topological sense and does maintain homoclinic/heteroclinic solutions as diffusive coupling coefficients tend to zero. Computational experiments are used to illustrate the dynamics of such solutions in the hepatitis C evolution model.

Keywords: hepatitis C model; homoclinic/heteroclinic solution; generalized differential operator; inverse balancing

MSC 2010: 34A34; 35C07; 92B05

1 Introduction

Modeling of biomedical processes using differential equations has become more and more widespread over recent years [6, 12, 29]. Various differential equation models on the use of oncolytic viruses as therapeutic agents against cancer are discussed in [28]. A clinically validated model of tumor-immune cell interactions is considered in [4]. A new mathematical model for the explanation of the failure of cancer chemotherapy treatment is presented in [22]. A mathematical model based on differential equations is used to describe the interactions between Ebola virus and wild-type Vero cells in vitro in [21].

Beginning with the classical paper by Neumann et al [20], various differential equation models for the modeling of hepatitis virus infection have been proposed. Global dynamics of a delay differential model of hepatitis B infection evolution are studied in [5, 27]. The transmission of hepatitis C virus (HCV) among injecting drug users is modeled using ordinary differential equations in [11]. A mathematical multi-scale

model of the within-host dynamics of HCV infection is used to study patients under treatment with direct acting antiviral medication in [3]. The authors of [2] give a review of recent HCV kinetics models.

Reluga et al [25] present the following model of hepatitis C virus infection that explicitly includes proliferation of infected and uninfected hepatocytes:

(1)Tt^=s^+rTT1−T+ITmax−dTT−1−ηβVT+q^I;It^=rII1−T+ITmax+1−ηβVT−dII−q^I;Vt^=1−ϵpI−cV,

where pt is time; Tt^represents uninfected hepatocytes; It^represents infected cells and Vt^represents free virus population. The parameters of (1) have the following meaning: β is the rate of infection per free virus per hepatocyte; c is the immune virus clearance rate; p is the free virus production rate per infected cell; dT, dI are death rates for uninfected hepatocytes and infected cells respectively; rT, rI are parameters of the logistic proliferation of T and I respectively; logistic proliferation happens only if T < T_max; parameters s^andq^represent the increase rate of uninfected hepatocytes through immigration and spontaneous cure by noncytolytic process respectively; finally the effect of antiviral treatment reduces the infection rate by a fraction η and the viral production rate by a fraction ϵ. Ranges of parameters are given in [25].

As shown in [25], V can be solved explicitly for patients in a steady state before treatment. Furthermore, introducing dimensionless state variables and parameters transforms (1) into:

(2)xτ=x1−x−y−1−θbxy+qy+s;yτ=ry1−x−y+1−θbxy−dy−qy,

where x, y are dimensionless state variables for uninfected hepatocytes and infected cells respectively; r, b, θ, d, q, s P R are real parameters.

System (2) can be rewritten in a more general form:

(3)xτ=a0+a1x+a2x2+a3xy+a4y;xτ=c=u;yτ=b0+b1y+b2y2+b3xy+b4x;xτ=c=v,

where c, u, v, a_k , b_k P R, k = 1, . . . , 4.

The main objective of this paper is to study soliton-like dynamics of the system (3). Note that since (3) is not a system of nonlinear partial differential equations (PDEs), soliton (or solitary) solutions cannot exist, due to their definition being closely connected to concrete physical phenomena. However, as is demonstrated in the paper, solutions that exhibit analogous dynamics to those observed in solitary solutions, can be constructed for system (3). Since the phase trajectories of these solutions are homoclinic or heteroclinic, we refer to such solutions and homoclinic/heteroclinic solutions.

In the case a₄ = b₄ = 0, system (3) has already been shown to admit homoclinic/heteroclinic solutions [19], [15]. Solutions described in [19] have simple monotonous transitions from two steady states, while those found in [15] exhibit much more complicated transient effects. Because of this reason, only the latter homoclinic and heteroclinic solutions to (2), (3) are considered.

Using the inverse balancing and generalized differential operator techniques, explicit homoclinic and heteroclinic solution existence conditions are obtained in terms of the parameters of (2). These conditions, together with explicit expressions of such solutions, provide insight not only into HCV model (2), but also other models of nonlinear evolution.

Note that the application of direct techniques to compute the homoclinic/heteroclinic trajectories of (3) is not straightforward. For example, computation of the first integral requires the solution of the following first-order ODE:

(4)yx=a0+a1x+a2x2+a3xy+a4yb0+b1y+b2y2+b3xy+b4x.

While the above ODE can be integrated for some parameter values, there is no general method to determine such cases. Furthermore, the generalized differential operator technique yields not only phase trajectories of (3), but also its general solution and the conditions with respect to a₀, . . . , a₄; b₀, . . . , b₄ under which homoclinic/heteroclinic solutions exist.

2 Preliminaries

2.1 Power series and their extensions

In this paper, functions of the following power series form are considered:

(5)f(z)=∑j=0+∞ajzjj!,

where z, a_j P C. The coefficients of power series (33) are constructed via generalized differential operator technique, described in the following sections of the paper.

We treat the convergence of series (33) as follows. If (33) converges in some ball |z| < R; R ¡ 0, then it is possible to extend (33) to a wider complex domain (not including the singularities of (33)) via classical extension techniques. Let t P R denote a real argument of this extended function. Inserting t into the extension of (33) yields a real power series f (x) defined for values not necessarily in the radius |t| < R. For the purposes of this paper, we consider f (x) and its power series representation to be congruent.

2.2 Monotonous and non-monotonous homoclinic/heteroclinic solutions

First, let us consider monotonous homoclinic and heteroclinic solutions of the following soliton-like form [23, 26]:

(6)xτ;c,u,v=σexpητ−c−x1expητ−c−τ1(x);

(7)yτ;c,u,v=γexpητ−c−y1expητ−c−τ1(y),

where η ≠ 0, σ, γ ∈ R are constants; τ1(x),τ1(y)depend on initial conditions u, v.

The biological interpretation of (6), (7) represents the transition from the size of population of cells before therapy to the size of population after therapy. However, this transition is monotonous; the solutions shown in Fig. 1 (a) describe the difference between the sizes of populations before and after therapy, and the transition between the steady states.

$Figure 1 Monotonous (a) and non-monotonous homoclinic/heteroclinic solutions (b). Black and gray lines represent x(τ) and y(τ) respectively. The parameters of solutions read: (η=σ=1;γ=1/2;c=0;x1=3;x2=1;τ1(x)=−2;τ2(x)=−5;y1=−8;y2=−1/2;τ1(y)=−3;τ2(y)=−1$(\eta = \sigma = 1; \gamma = 1/2; c = 0; x_1 = 3; x_2 = 1; \tau^{(x)}_1 = -2; \tau^{(x)}_2 = -5; y_1 = -8; y_2 = -1/2; \tau^{(y)}_1 = -3; \tau^{(y)}_2 = -1$.$

Figure 1

Monotonous (a) and non-monotonous homoclinic/heteroclinic solutions (b). Black and gray lines represent x(τ) and y(τ) respectively. The parameters of solutions read: (η=σ=1;γ=1/2;c=0;x1=3;x2=1;τ1(x)=−2;τ2(x)=−5;y1=−8;y2=−1/2;τ1(y)=−3;τ2(y)=−1.

Non-monotonous homoclinic/heteroclinic solutions read [7, 26]:

(8)xτ;c,u,v=σexpητ−c−x1expητ−c−x2expητ−c−τ1(x)expητ−c−τ2(x);

(9)yτ;c,u,v=γexpητ−c−y1expητ−c−y2expητ−c−τ1(y)expητ−c−τ2(y),

where η ≠ 0, σ, γ ∈ R are constants; τk(x),τk(y),xk,yk,k=1,2depend on initial conditions u, v.

Solutions (8), (9) describe much more complex transition processes between the steady states. The size of the population of cells during the transient process exceeds populations both at the beginning and the end of the therapy if only the considered solutions have minimum points (the black line in Fig. 1 (b)) Analogously,

solutions with maximum points describe complex transitions from the population of cells before and after the treatment (the gray line in Fig. 1 (b))

From the biological point of view, transient processes governed by homoclinic and heteroclinic solutions highlight important phenomena. Let us consider the dynamics of uninfected cells (the black line in Fig. 1 (b)) The population of uninfected cells after the therapy becomes lower than the population before the therapy. However, the number of uninfected cells grows during the therapy and exceeds the population of uninfected cells at the beginning of the computational experiment (Fig. 1 (b))

Note that the negative values of cell population x(τ) and y(τ) are a consequence of the non-dimensionalization of system (1).

2.3 Solution transformation

In the following derivations, the standard independent variable transformation will be used:

(10)t:=expητ−c;c^:=expητ−c.

Using (10), homoclinic/heteroclinic solutions (8), (9) can be written as:

(11)x^t;c,u,v:=x1ηln⁡τ;1ηln⁡c,u,v=σt−x^1t−x^2t−t1(x)t−t2(x);

(12)y^τ;c,u,v:=y1ηln⁡τ;1ηln⁡c,u,v=γt−y^1t−y^2t−t1(y)t−t2(y),

where tk(x)=c^τk(x),tk(y)=c^τk(y),x^k=c^xk,y^k=c^yk,k=1,2Using partial fractions (11), (12) can be rewritten as:

(13)x^=σ+λ11−ρ1t−c^+λ21−ρ2t−c^;

(14)y^=γ+μ11−ν1t−c^+μ21−ν2t−c^,

where λ_k , μ_k , ρ_k , ν_k , k = 1, 2 are functions of u, v.

2.4 Generalized differential operator technique

In this section, a summary on the generalized differential operator technique for the construction of solutions to ordinary differential equations in presented. More detailed derivations can be found in [16].

2.4.1 Generalized differential operators

Let P(c, u, v), Q(c, u, v) be trivariate analytic functions. A generalized differential operator D_cuv reads:

(15)Dcuv:=Dc+P(c,u,v)Du+Q(c,u,v)Dv,

where Dβ:=∂∂βfor any variable β. Standard properties of differentiation operators hold true for (15) [14]:

(16)Dcuvfgc,u,v=dfdgDc,uvg;

(17)Dcuvfg=gDcuvf+fDcuvg;

(18)Dcuvfg=gDcuvf−fDcuvgg2,

where f , g denote arbitrary functions analytic in c, u, v.

2.4.2 Multiplicative operators

Using (15), the multiplicative operator can be constructed:

(19)M:=∑j=0+∞tjj!Dcuvj,

where t is an arbitrary real variable. Operator (19) has two important properties:

(20)Mcm=t+cm,m=0,1,…;

(21)Mfc,s,t=ft+c,Ms,Mt.

Note that (20) follows immediately from the definition of (19). Without loss of generality, the proof of (21) for multiplicative operator M=∑j=0+∞tjj!P(u,v)Du+Q(u,v)Dvjis presented below.

Let y₁ :═ Mu ═ y₁(t, u, v), y₂ :═ Mv ═ y₂(t, u, v), z :═ Mf (u, v) ═ z(t, u, v) and w :═ f (Mu,Mv) ═ f (y₁, y₂). To prove (21), it needs to be shown that z ═ w for all t, u, v.

Note that:

(22)Dtz=Dt∑j=0+∞tjj!Duvjf(u,v)=∑j=1+∞tj−1(j−1)!Duvjf(u,v)=DuvMf(u,v)=Duvz=PDuz+QDvz.

Thus, the function z(t, u, v) satisfies the partial differential equation:

(23)∂z∂t=P∂z∂u+Q∂z∂v,

with initial condition z(0, u, v) ═ f (u, v) that follows from the definition of z.

Analogously, it is shown that:

(24)∂yk∂t=P∂yk∂u+Q∂yk∂v; k=1,2;

with y₁(0, u, v) ═ u and y₂(0, u, v) ═ v. Using (24) and the definition of w yields:

Note that w satisfies the initial condition w(0, u, v) ═ f (y₁(0, u, v), y₂(0, u, v)) ═ f (u, v), thus z and w coincide, which results in the proof of (21).

Construction of general solutions to ODEs requires one final operator which is denoted as the generalized multiplicative operator:

(26)G:=∑j=0+∞(t−c)jj!Dcuvj.

Operator G has two properties analogous to (20), (21):

(27)Gcm=tm,m=0,1,…;

(28)Gf(c,u,v)=fx,Gu,Gv,

where f is a trivariate analytic function. The proof of (28) follows from (21):

(29)Mf(c,u,v)=ft+c,∑j=0+∞tjj!Dcuvju,∑j=0+∞tjj!Dcuvjv.

Substituting t for t ═ c yields (28).

2.4.3 Construction of solutions to ODEs

Let us consider the following system of ODEs:

(30)x^t=P(t,x^,y^); x^|t=c^=u;y^t=Q(t,x^,y^); y^|t=c^=v,

where P, Q are analytic functions. The generalized differential operator respective to (30) reads [13]:

(31)Dc^uv:=Dc^+Pc^,u,vDu+Qc^,u,vDv.

Using (31), general solution to (30) is expressed as [13, 14]:

(32)x^=Gu=∑j=0+∞t−c^jj!Dc^uvju;y^=Gv=∑j=0+∞t−c^jj!Dc^uvjv.

The convention Dc^uv0=Iwhere I is the identity operator, is used.

Identities (32) can be proven using properties (21) and (28) derived in the previous section. Consider operators M, G defined with respect to the generalized differential operator (31). First, let z=z(t,c^,u,v)=Mu and w=w(t,c^,u,v)=MvProperty (21) yields:

(33)Dtz=DtMu=Dc^uvMu=MDc^uv=MPc^,u,v=Pt+c,Mu,Mv=P(t+c,z,w).

Analogously,

(34)Dtw=Q(t+c,z,w).

Selecting x^=z(t−c,c^,u,v)=Guandy^=w(t−c,c^,u,v)=Gvyields the system (30). Furthermore, the definition of operator G yields that x^|t=c^=u and y^ |t=c^=v,thus (32) hold true.

In the following derivations,

(35)pj=pjc^,u,v:=Dc^uvju;qj=qjc^,u,v:=Dc^uvjv;j=0,1,…

will be used, which transforms (32) into:

(36)x^=∑j=0+∞t−c^jj!pj;

(37)y^=∑j=0+∞t−c^jj!qj.

Furthermore, coefficients p _j , q_j satisfy recurrence relations:

(38)pj+1=Dc^uvpj;qj+1=Dc^uvqj.

3 Existence of homoclinic/heteroclinic solutions in (30)

Let ρ₁ ≠ ρ₂. If (30) admits solutions (13), (14) then (13) and (36) must be equal. Expanding (13) in a power series and equating to (36) yields:

(39)σ+λ1+λ2+∑j=1+∞t−c^jj!j!λ1ρ1j+j!λ2ρ2j=∑j=0+∞t−c^jj!pj.

Note that p₀ ═ u by (35), thus (39) yields:

(40)p0=u;

(41)pj=j!λ1ρ1j+λ2ρ2j,j=1,2,….

Analogous derivations with respect to y and ν₁ ≠ ν₂ result in:

(42)q0=v;

(43)qj=j!μ1ν1j+μ2ν2j,j=1,2,….

Thus (30) admits solutions (13), (14) if and only if (41), (43) hold true.

Theorem 3.1System (30) admits homoclinic/heteroclinic solutions (13), (14) with ρ₁ ≠ ρ₂if and only if

(44)λk=p2−2ρlp12ρkρk−ρl;μk=q2−2νlq12ρkνk−νl;

(45)Dc^uvρk=ρk2;Dc^uvνk=νk2;

(46)Dc^uvλk=λkρk;Dc^uvμk=μkνk;

(47)3p12p42−36p1p2p3p4+32p1p33+36p23p4−36p22p32≠0;

(48)3q12q42−36q1q2q3q4+32q1q33+36q23q4−36q22q32≠0,

k, l = 1, 2; k ≠ l.

Proof. It will be proven that (41), (43) hold true if and only if (44)–(48) hold true.

Necessity Let (41) hold true. Taking j = 1, 2 yields:

(49)p1=λ1ρ1+λ2ρ2;

(50)p2=2λ1ρ12+λ2ρ22.

Solving the above equations for λ₁, λ₂ results in (44).

Equation (41) yields the following determinant equality:

(51)detp11!p22!p33!p22!p33!p44!1ρkρk2=0;k=1,2.

Expanding the left side of (51) yields:

(52)Δ2ρk2−Δ1ρk+Δ0=0;k=1,2,

where

(53)Δ2=p1p33!−p22!2;

(54)Δ1=p1p44!−p2p32!⋅3!;

(55)Δ0=p2p42!⋅4!−p33!2.

Solving (52) for ρk results in: _k

(56)ρ1,2=Δ1±Δ12−4Δ2Δ02Δ2.

Since ρ₁ ≠ ρ₂, the discriminant Δ12−4Δ2Δ0≠0which results in condition (47).

Denoting Θ:=Δ12−4Δ2Δ0and applying operator Dc^uvto (56) results in:

(57)Dc^uvρ1,2=±Θ−Δ1Dc^uvΔ2+Δ2Dc^uvΔ1±Dc^uvΘ2Δ22.

Using recursion (38) it can be obtained that:

(58)Dc^uvΔ2=p1p46−p2p33;

(59)Dc^uvΔ1=p1p524−p2p424−p3212;

(60)Dc^uvΔ0=p2p548−5p3p4144;

(61)Dc^uvΘ=1ΘΔ1Dc^uvΔ1−2Δ0Dc^uvΔ2−2Δ2Dc^uvΔ0.

Relation (41) transforms (53)–(55) and 𝛩 into:

(62)Δ2=λ1λ2ρ1ρ2ρ1−ρ22;

(63)Δ1=λ1λ2ρ1ρ2ρ1+ρ2ρ1−ρ22;

(64)Δ0=5λ1λ2ρ12ρ22ρ1−ρ22;

(65)Θ=λ1λ2ρ1ρ2ρ1−ρ23.

Furthermore,

(66)Dc^uvΔ2=4λ1λ2ρ1ρ2ρ1+ρ2ρ1−ρ22;

(67)Dc^uvΔ1=λ1λ2ρ1ρ2ρ1−ρ225ρ12+8ρ1ρ2+5ρ22;

(68)Dc^uvΔ0=5λ1λ2ρ12ρ22ρ1+ρ2ρ1−ρ22;

(69)Dc^uvΘ=5λ1λ2ρ1ρ2ρ1+ρ2ρ1−ρ23.

Inserting (62)–(69) into (57) yields (45).

Applying operator Dc^uvto (44) and using (45) yields:

(70)Dc^uvλk=4p1ρkρl−p22ρk+3ρl+p32ρkρk−ρl.

Inserting (41) into (70) results in (46).

Sufficiency Condition (44) yields:

(71)p1=λ1ρ1+λ2ρ2.

Applying operator Dc^uvto (71) results in:

(72)Dc^up1=ρ1Dc^uvλ1+λ1Dc^uvρ1+ρ2Dc^uvλ2+λ2Dc^uvρ2=2λ1ρ12+λ2ρ22=p2.

Continuing by induction yields (41).

The proof for parameters of y is analogous.

Corollary 3.1If conditions of Theorem 3.1 hold true, then the third and higher order Hankel determinants of sequencespjj!,qjj!;j=1,2,…are equal to zero

(73)Hp(n)=detpj+k−2(j+k−2)!1≤j,k≤n+1=0;

(74)Hq(n)=detqj+k−2(j+k−2)!1≤j,k≤n+1=0,

n = 3, 4, . . .

Proof. Proof results from the derivation of Theorem 3.1 and (41), (43).

4 Necessary homoclinic/heteroclinic solution existence conditions in (3)

The inverse balancing technique can be used to determine necessary existence conditions of solutions (8), (9) to (3). The main principle of this technique is to insert the solution ansatz into the considered equations and obtain a system of equations linear in system parameters a_k , b _k , k = 0, . . . , 4. The inverse balancing technique has been successfully used to obtain necessary solution existence conditions in a variety of nonlinear ordinary and partial differential equations [10, 15, 18]. Note that the inverse balancing technique does not possess the drawbacks associated with various solution construction (or direct ansatz) methods, which have attracted a significant amount of criticism [1, 8, 9, 17, 24].

4.1 Transformation of (3)

Using the substitution (10), system (3) is transformed to:

(75)ηtx^t=a0+a1x^+a2x^2+a3x^y^+a4y^;ηty^t=b0+b1y^+b2y^2+b3x^y^+b4x^,

with initial conditions

(76)x^|t=c^=u; y^|t=c^=v.

The following notations are introduced:

(77)X(t):=t−x^1t−x^2;Y(t):=t−y^1t−y^2;

(78)Tx(t):=t−t1(x)t−t2(x);Ty(t):=t−t1(y)t−t2(y),

which transform solutions (11), (12) to:

(79)x^=σX(t)Tx(t);y^=γY(t)Ty(t).

4.2 Necessary existence conditions for (79) in (75)

Following the inverse balancing technique, solution ansatz (79) is inserted into (75). After simplification, (75) reads:

(80)ηtσTyXtTx−XTxt=a0Tx2Ty+a1σTxTyX+a2σ2X2Ty+a3σγXYTx+a4γYTx2;

(81)ηtγTxYtTy−YTyt=b0Ty2Tx+b1γTyTxY+b2γ2Y2Tx+b3σγYXTy+b4γXTy2.

Equation (78) results in:

(82)Txt1(x)=Txt2(x)=Tyt1(y)=Tyt2(y)=0;

(83)(Tx)t|t=t1(x)=t1(x)−t2(x); (Tx)t|t=t2(x)=t2(x)−t1(x)

(84)(Ty)t|t=t1(y)=t1(y)−t2(y); (Ty)t|t=t2(y)=t2(y)−t1(y).

Letting t=t1(x),t2(x)in (80), t=t1(y),t2(y)in (81) and using (82)–(84) yields the following equations:

(85)Tyt1(x)a2σ2X2t1(x)+ησt1(x)t1(x)−t2(x)Xt1(x)=0;

(86)Tyt2(x)a2σ2X2t2(x)+ησt2(x)t2(x)−t1(x)Xt2(x)=0;

(87)Txt1(y)b2γ2Y2t1(y)+ηγt1(y)t1(y)−t2(y)Yt1(y)=0;

(88)Txt2(y)b2γ2Y2t2(y)+ηγt2(y)t2(y)−t1(y)Yt2(y)=0.

Equations (85)–(88) have nontrivial solutions only if:

(89)t1(x)=t1(y);t2(x)=t2(y),

thus (75) (and conversely (3)) only admits homoclinic/heteroclinic solutions with equal denominators. Let t1:=t1(x)=t1(y)andt2:=t2(x)=t2(y) Equation (89) transforms (79) into:

(90)x^=σX(t)T(t);y^=γY(t)T(t),

where T(t) := (t = (₁) (t = t₂).

4.3 Necessary existence conditions for (90) in (75)

If (89) holds true, (80), (81) read:

(91)ηtσXtT−XTt=a0T2+a1σXT+a2σ2X2+a3σγXY+a4γYT;

(92)ηtγYtT−YTt=b0T2+b1γYT+b2γ2Y2+b3σγXY+b4σXT.

Note that

(93)Tt1=Tt2=Xx^1=Xx^2=Yy^1=Yy^2=0,

and

(94)Tt=2t−t1−t2;Xt=2t−x^1−x^2;Yt=2t−y^1−y^2.

Taking t = t₁, t₂ in (91) and using (93), (94) yields:

(95)ηt1t2−t1=σXt1a2+γYt1a3;

(96)ηt2t1−t2=σXt2a2+γYt2a3.

Analogous computations with respect to (92) result in:

(97)ηt1t2−t1=γYt1b2+σXt1b3;

(98)ηt2t1−t2=γYt2b2+σXt2b3.

Solution of (95)–(98) with respect to a₂, a₃, b₂, b₃ reads:

(99)a2=b3=ησt2−t1t1Yt2+t2Yt1Xt1Yt2−Xt2Yt1;

(100)b2=a3=ηγt1−t2t2Xt1+t1Xt2Xt1Yt2−Xt2Yt1.

Similarly, taking t=x^1,x^2in (91) and t=y^1,y^2in (92) yields the following solutions for a₀, a₄, b₀, b₄:

(101)a0=ησx^1−x^2x^1Yx^2+x^2Yx^1Tx^1Yx^2−Tx^2Yx^1;

(102)a4=ησx^2−x^1x^1Tx^2+x^2Tx^1γTx^1Yx^2−Tx^2Yx^1;

(103)b0=ηγy^1−y^2y^1Xy^2+y^2Xy^1Ty^1Xy^2−Ty^2Xy^1;

(104)b4=ηγy^2−y^1y^1Ty^2+y^2Ty^1σTy^1Xy^2−Ty^2Xy^1.

Finally, taking t = 0 in (91), (92) yields a₁, b₁:

(105)a1=−1σa0+σ2a2+σγa3+γa4;

(106)b1=−1γb0+γ2b2+σγb3+σb4.

Note that there are 10 parameters in (75) and (91), (92) yields a non-degenerate system of 10 linear balancing equations, thus no constraints on the parameters of solution (90) needs to be imposed. However, as shown by (99), (100) conditions a₃ = b₂ and b₃ = a₂ must hold if (75) admits solution (90).

The results of this section are summarized in the following Lemma.

Lemma 4.1System (3) admits homoclinic/heteroclinic solutions (8), (9) only if

(107)τ1(x)=τ1(y);τ2(x)=τ2(y);

(108)a3=b2;b3=a2.

Note that condition (107) results from (89) and substitution (10). Also, ρ_k = ν_k; k = 1, 2 in (13), (14) when (107) holds true.

5 Construction of homoclinic/heteroclinic solutions to (3)

In this section, explicit expressions of homoclinic and heteroclinic solutions to (3) are constructed. It is assumed that the necessary existence conditions (107), (108) hold true.

5.1 Derivation of parameter η

Parameter η is derived using Corollary 3.1. Consider the following Hankel determinants:

(109)Hp(3)=p11!p22!p33!p22!p33!p44!p33!p44!p55!;Hq(3)=q11!q22!q33!q22!q33!q44!q33!q44!q55!.

Parameter η must be chosen to satisfy

(110)Hp(3)=0;Hq(3)=0.

Furthermore, η can only depend on coefficients a₀, . . . , a₄; b₀, . . . , b₄, otherwise Theorem 3.1 does not hold true and obtained solutions would not be valid for all initial conditions.

It can be observed that:

(111)Hp(3)=1η9c^9A6(u,v)η6+A4(u,v)η4+A2(u,v)η2+A0(u,v);

(112)Hq(3)=1η9c^9B6(u,v)η6+B4(u,v)η4+B2(u,v)η2+B0(u,v).

Thus, roots of equations (111), (112) with respect to η that do not depend on u, v must be found. Note that:

(113)A6=121601/3K3,K:=u2a2+a3v+a1u+a4v+a0;

and

(114)A4=Fu,vK,

where F is a polynomial in u, v.

Since the roots η must not depend on initial conditions, any values of u, v can be chosen and inserted into (111). Let

(115)v=f(u)=−a2u2+a1u+a0b2u+a4,

then A₆ = A₄ = 0 and using (111), η² can be expressed as:

(116)η2=−A0(u,f(u))A2(u,f(u)).

The numerator and denominator of (116) depend linearly on u:

(117)η2=α1u+α0β1u+β0,

where α_k , β_k are functions of a₀, . . . , a₄; b₀, . . . , b₄.

Analogous computations with respect to Hq(3)lead to:

(118)u=g(v)=−b2v2+b1v+b0a2v+b4;

(119)η2=−B0(g(v),v)B2(g(v),v)=α^1v+α^0β^1v+β^0.

Parameter η does not depend on u, v only if:

(120)α1β0−α0β1=0;

(121)α^1β^0−α^0β^1=0

Note that:

(122)1a0b22−a1a4b2+a2a42α1β0−α0β1=1b0a22−b1b4a2+b2b42α^1β^0−α^0β^1,

which leads to the following sufficient existence condition for homoclinic/heteroclinic solutions to (3):

(123)9a0a1a2+9b0b1b2−18a0a2b1−18b0b2a1+3a1b12+3b1a12−2a13−2b13−9a1a4b4−9b1b4a4+27a0b2b4+27b0a2a4=0.

If (123) holds true, η can be computed from either (117) or (119). Furthermore, if (123) holds true, the parameter η does not depend on initial conditions c, u, v.

5.2 Necessary and sufficient existence conditions for homoclinic/heteroclinic solutions to (3)

Theorem 3.1, Lemma 4.1 and condition (123) together with computer algebra computations result in the following theorem.

Theorem 5.1System (3) admits homoclinic/heteroclinic solutions (8), (9) if and only if conditions (107), (108) and (123) hold true.

Note that

(124)ρk=ρk∗c^;k=1,2,

where ρk∗=ρk∗u,v

Relations between parameters of (13), (14) and (8), (9) read:

(125)τk=1+1ρk∗;k=1,2;

(126)x1,2=12Ax±Ax2−4Bx;

(127)y1,2=12Ay±Ay2−4By,

where

(128)Ax:=λ1σρ1∗+λ2σρ2∗+τ1+τ2;

(129)Bx:=τ1τ2+λ1τ2σρ1∗+λ2τ1σρ2∗;

(130)Ay:=μ1γρ1∗+μ2γρ2∗+τ1+τ2;

(131)By:=τ1τ2+μ1τ2γρ1∗+μ2τ1γρ2∗.

Parameters σ,γ read:

(132)σ=u−λ1−λ2;γ=v−μ1−μ2.

Note that (117) yields two values for η, however, it is sufficient to consider only the positive or negative root of (117) to obtain the general solution to (3) when Theorem 5.1 holds true, because the sign of η can be interchanged:

(133)x=σexpητ−c−x1expητ−c−x2expητ−c−τ1expητ−c−τ2=σx1x2exp−ητ−c−1x1exp−ητ−c−1x2τ1τ2exp−ητ−c−1τ1exp−ητ−c−1τ2.

As demonstrated in [15], the value x1x2τ1τ2does not depend on initial conditions, which proves that changing the sign of η does not yield new solutions.

6 Homoclinic/heteroclinic solutions to hepatitis C model (2)

6.1 Existence conditions

Comparing (2) to (3) it can be observed that:

(134)a0=s;a1=1;a2=−1;a3=bθ−1−1;a4=q;

(135)b0=0;b1=r−d−q;b2=−r;b3=b1−θ−r;b4=0.

To preserve biological significance of system (2), the parameters (134), (135) must satisfy q, s, r ≥ 0; b P [0^-2; 10³]; d ∈ [10^-3; 10²][25].

Using Theorem 5.1 conditions for the existence of homoclinic/heteroclinic solutions to (2) can be derived. Note that only homoclinic/heteroclinic solutions with τk(x)=τk(y)can be considered. Inserting (134), (135) into (108) yields two congruent equations:

(136)1−θb−r=−1;bθ−1−1=−r.

Both equations are satisfied if parameter r reads:

(137)r=b1−θ+1.

Let (137) hold true. Denote := b (1 ‒ θ) + 1 ‒ d ‒ q = r ‒ d ‒ q. Inserting (134), (135) into condition (123) yields:

(138)s=19Ω2−Ω−2.

Equations (137) and (138) result in the following corollary.

Corollary 6.1Hepatitis C model (2) admits homoclinic/heteroclinic solutions if and only ifτk(x)=τk(y)(137) and (138) hold true.

Computer algebra computations prove that when Corollary (6.1) holds true, parameters y₁ = y₂ = 0.

6.2 Equilibria

Let (137) and (138) hold true. The equilibria of (2) read:

(139)x∗(1)=23−Ω3;y∗(1)=0;

(140)x∗(2)=13+Ω3;y∗(1)=0;

(141)x∗(3)=1−θs+qb−q2+1−dq+s1−θ2b2−1−θd+q−1b−d;

(142)y∗(3)=22Ω−1291−θ2b2−1−θd+q−1b−d.

Equilibrium point (141), (142) is a stable node as τ ➝ +∞:

(143)limτ→+∞xτ,yτ=σ,γ.

Equilibrium point (139) is an unstable node as τ ➝ ‒∞:

(144)limτ→−∞xτ,yτ=σx1x2τ1τ2,0.

The remaining equilibrium point (140) is a saddle point.

6.3 Computational experiment

Let us consider the following system:

(145)xτ=x(1−x−y)+18xy+4y+49; x|τ=c=u;

(146)yτ=19y(1−x−y)−18xy−16y; y|τ=c=v.

The above system corresponds to (2) with the following parameters:

(147)b=24;θ=14;q=4;d=12;r=19;s=49.

Note that parameters (147) satisfy the guidelines given in [25] for biologically significant systems. Further more, conditions of Corollary 6.1 are satisfied, thus homoclinic/heteroclinic solutions to (145), (146) do exist.

Equation (117) yields:

(148)η=±53.

As noted previously, it is sufficient to consider one value of η to obtain the general solution to (145), (146). In subsequent computations the value η=53is used.

Theorem 3.1 yields the following parameters of homoclinic/heteroclinic solutions:

(149)ρ1,2=110c^−3u−57v−1±Φ;

(150)λ1,2=3Φ33−1083v2+57u−202v+5u+53Φ∓3648v2∓570u−760v∓u−43ΦΦ9u+171v+3∓3Φ;

(151)μ1,2=v3u+57v−9Φ∓A−30u−192v+40Φ3u+57v+1−Φ,

$Figure 2 Homoclinic/heteroclinic solutions to (145), (146). Black and gray lines correspond to x and y respectively. Dotted lines denote singularity points in (c) and (d). Initial conditions are u = ‒1, v = 1{10 in (a); u = 105{100, v = 3{100 in (b); u = ‒2, v = ‒1{100 in (c); u = ‒3; v ‒ 2{100 in (d). Labels (a), (b), (c), (d) correspond to respectively labeled phase plane trajectories in Fig. 3.$

Figure 2

Homoclinic/heteroclinic solutions to (145), (146). Black and gray lines correspond to x and y respectively. Dotted lines denote singularity points in (c) and (d). Initial conditions are u = ‒1, v = 1{10 in (a); u = 105{100, v = 3{100 in (b); u = ‒2, v = ‒1{100 in (c); u = ‒3; v ‒ 2{100 in (d). Labels (a), (b), (c), (d) correspond to respectively labeled phase plane trajectories in Fig. 3.

$Figure 3 Phase portrait of (145), (146). Gray circles denote the stable and unstable nodes; diamond denotes the saddle point. Solid black lines correspond to solution trajectories. Dashed gray parabola corresponds to the separatrix between solutions with elliptic and hyperbolic trajectories. Dashed gray lines denote stable and unstable manifolds of the saddle point. Labels (a), (b), (c), (d) correspond to respective parts of Fig. 2. Trajectories in the solid gray and horizontally striped filled regions are elliptic and hyperbolic respectfully and do not have singularities. Trajectories in the unfilled regions are hyperbolic and have one singularity. Trajectories in vertically striped regions are hyperbolic and have two singularities.$

Figure 3

Phase portrait of (145), (146). Gray circles denote the stable and unstable nodes; diamond denotes the saddle point. Solid black lines correspond to solution trajectories. Dashed gray parabola corresponds to the separatrix between solutions with elliptic and hyperbolic trajectories. Dashed gray lines denote stable and unstable manifolds of the saddle point. Labels (a), (b), (c), (d) correspond to respective parts of Fig. 2. Trajectories in the solid gray and horizontally striped filled regions are elliptic and hyperbolic respectfully and do not have singularities. Trajectories in the unfilled regions are hyperbolic and have one singularity. Trajectories in vertically striped regions are hyperbolic and have two singularities.

where

(152)Φ=Φ(u,v):=9u2+342uv+3249v2+6u−642v+1.

Derivations given in Subsection 5.2 result in:

(153)σ=92189;γ=29189;

(154)τ1,2=Φ∓3u+57v−9Φ∓3u+57v+1;

(155)x1,2=v3138−15u+93v−5±Φ;y1,2=0.

Solutions with parameters (153)–(155) are pictured in Fig. 2. Note that there are three types of solutions – non-singular solutions (a), (b); solutions with one singularity (c) and solutions with two singularities (d).

The phase plane of (145), (146) can be seen in Fig. 3. Note that labels (a), (b), (c), (d) on the phase plane correspond to respectively labeled solutions pictured in Fig. 2. System(145), (146) has the following equilibria 92189,29189- stable node; −13,0- unstable node; 43,0- saddle point.

It has been proven in [15] that homoclinic/heteroclinic solutions of the form (8), (9) correspond to phase plane trajectories that satisfy the general conic section equation:

(156)Ax2+Bxy+Cy2+Ex+Fy+G=0;A,B,C,E,F,G∈R.

Solution Fig. 3 (a) corresponds to an elliptic trajectory, while the remaining (b), (c), (d) have hyperbolic trajectories. Furthermore, there is a single solution that satisfies the parabola equation:

(157)Φ(x,y)=0.

Curve (157) is a separatrix that separates solutions with and without singularities in the phase plane (see dashed gray parabola in Fig. 3).

Stable and unstable manifolds of the saddle point are obtained by setting the numerator and denominator of τ_1,2 to zero [15]. This yields that the stable manifold of the saddle point is the x-axis, while the unstable manifold lies on the straight line y=−532x+524.Manifolds of the saddle point correspond to dashed gray lines in Fig. 3.

7 Concluding remarks

Homoclinic and heteroclinic solutions to hepatitis C evolution model (2) have been constructed in this paper. Inverse balancing and generalized differential operator techniques have enabled the derivation of explicit necessary and sufficient homoclinic and heteroclinic solution existence conditions with respect to the parameters of system (2). Furthermore, it has been shown that these existence conditions are satisfied when (2) described a biologically significant system of HCV evolution.

It has been demonstrated that transient processes of the derived solutions to (2) reveal important phenomena for understanding hepatitis C virus infection dynamics. Even though antiviral therapy reduces the number of infected cells (comparing the beginning to the end of treatment), due to the transient processes during the therapy, population size of infected cells is higher than before or after therapy – if only the considered solutions are heteroclinic with maxima. Analogous biological interpretations can be made for heteroclinic solutions with minima. The population of healthy cells is lower than before or after treatment during antiviral therapy – if the number of uninfected hepatocytes is described by a heteroclinic solution possessing minima.

The main mathematical advancements of this paper can be characterized by new applications of inverse balancing technique and the development of generalized differential operator method for the solution of coupled differential equations with multiplicative and diffusive terms. As noted in Section 4, direct balancing techniques may yield wrong solutions; inverse balancing of such a complex system of nonlinear differential equations poses a number of technical problems. On the other hand, derivation of closed-form homoclinic/heteroclinic solutions and explicit conditions of their existence poses serious mathematical challenges. One of the main contributions of this paper are the necessary and sufficient conditions for the existence of these solutions in the hepatitis C evolution model.

Comparing the results of this paper with [15] it can be concluded that system (3) (and, by extension (2)) is structurally stable in the topological sense – when a₄, b₄ tend to zero, the phase plane continuously converges to the phase plane described in [15]. Moreover, structural stability can also be observed in homo-clinic/heteroclinic solution existence condition (123) – in the case a₄, b₄ → 0, such solutions also exist and the condition (123) is maintained. Since such effects are observed in systems with biological significance, they provide valuable insight not only into (2) but also other nonlinear evolution models.

Acknowledgement

This research was funded by a grant (No. MIP078/2015) from the Research Council of Lithuania. This research was also funded by Jiangsu Provincial Recruitment Program of Foreign Experts (Type B, Grant 172 no. JSB2017007).

References

[1] Aslan Í, Marinakis V., Some remarks on Exp-function method and its applications, Commun Theor Phys 2011, 56, 397–403.10.1088/0253-6102/56/3/01Search in Google Scholar

[2] Chatterjee A., Guedj J., Perelson A.S., Mathematical modeling of HCV infection: what can it teach us in the era of direct antiviral agents? Antivir Ther 2012, 17, 1171–1182.10.3851/IMP2428Search in Google Scholar PubMed PubMed Central

[3] Clausznitzer D., Harnisch J., Kaderali L., Multi-scale model for hepatitis C viral load kinetics under treatment with direct acting antivirals, Virus Res 2016, 218, 96–101.10.1016/j.virusres.2015.09.011Search in Google Scholar PubMed

[4] de Pillis L.G., Radunskaya A.E., Wiseman C.L., A validated mathematical model of cell-mediated immune response to tumor growth, Cancer Res 2005, 65, 7950–7958.10.1158/0008-5472.CAN-05-0564Search in Google Scholar PubMed

[5] Gourley S.A., Kuang Y., Nagy J.D., Dynamics of a delay differential equation model of hepatitis B virus infection, J Biol Dyn 2:140–153, 2008.10.1080/17513750701769873Search in Google Scholar PubMed

[6] Hoppensteadt F.C., Peskin C., Eds., Modeling and Simulation in Medicine and the Life Sciences 2002, Springer, New York.10.1007/978-0-387-21571-6Search in Google Scholar

[7] Kudryashov N.A., On new travelling wave solutions of the KdV and the KdV-Burgers equation, Commun Nonlinear Sci Numer Simul 2009, 14, 1891–1900.10.1016/j.cnsns.2008.09.020Search in Google Scholar

[8] Kudryashov N.A., Seven common errors in finding exact solutions of nonlinear differential equations, Commun Nonlinear Sci Numer Simul 2009, 14, 3507–3529.10.1016/j.cnsns.2009.01.023Search in Google Scholar

[9] Kudryashov N.A., Loguinova N.B., Be careful with Exp-function method, Commun Nonlinear Sci Numer Simul 2009, 14, 1891–1900.10.1016/j.cnsns.2008.07.021Search in Google Scholar

[10] Marcinkevicius R., Navickas Z., Ragulskis M., Telksnys T., Solitary solutions to a relativistic two-body problem, Astrophys Space Sci 2016, 361, 201.10.1007/s10509-016-2792-2Search in Google Scholar

[11] Martin N.K., Pitcher A.B., Vickerman P., Vassall A., Hickman M., Optimal control of hepatitis C antiviral treatment programme delivery for prevention amongst a population of injecting drug users, Plos One 2011, 6, e22309.10.1371/journal.pone.0022309Search in Google Scholar PubMed PubMed Central

[12] Mierke C.T., Ed., Physics of Cancer 2015, IOP Publishing, Bristol.10.1088/978-0-7503-1134-2Search in Google Scholar

[13] Navickas Z., Bikulciene L., Expressions of solutions of ordinary differential equations by standard functions, Math Model Anal 2006, 11, 399–412.10.3846/13926292.2006.9637327Search in Google Scholar

[14] Navickas Z., Bikulciene L., Ragulskis M., Generalization of Exp-function and other standard function methods, Appl Math Comput 2010, 216, 2380–2393.10.1016/j.amc.2010.03.083Search in Google Scholar

[15] Navickas Z., Marcinkevicius R., Telksnys T., Ragulskis M., Existence of second order solitary solutions to Riccati differential equations coupled with a multiplicative term, IMA J Appl Math 2016, 81, 1163–1190.10.1093/imamat/hxw050Search in Google Scholar

[16] Navickas Z., Ragulskis M., How far can one go with the Exp-function method? Appl Math Comput 2009, 211, 522–530.10.1016/j.amc.2009.01.074Search in Google Scholar

[17] Navickas Z., Ragulskis M., Bikulciene L., Be careful with the Exp-function method – additional remarks, Commun Nonlinear Sci Numer Simul 2010, 15, 3874–3886.10.1016/j.cnsns.2010.01.032Search in Google Scholar

[18] Navickas Z., Ragulskis M., Telksnys T., Existence of solitary solutions in a class of nonlinear differential equations with polynomial nonlinearity, Appl Math Comput 2016, 283, 333–338.10.1016/j.amc.2016.02.049Search in Google Scholar

[19] Navickas Z., Vilkas R., Telksnys T., Ragulskis M., Direct and inverse relationships between Riccati systems coupled with multiplicative terms, J Biol Dyn 2016, 10, 297–313.10.1080/17513758.2016.1181801Search in Google Scholar PubMed

[20] Neumann A.U., Lam N.P., Dahari H., Gretch D.R.,Wiley T.E., Layden T.J., Perelson A.S., Hepatitis C viral dynamics in vivo and the antiviral efficacy of interferon-alpha therapy, Science 1998, 282, 103–107.10.1126/science.282.5386.103Search in Google Scholar PubMed

[21] Nguyen V.K., Binder S.C., Boianelli A., Meyer-Hermann M., Hernandez-Vargas E.A., Ebola virus infection modeling and identifiability problems, Front Microbiol 2015, 6, 257.10.3389/fmicb.2015.00257Search in Google Scholar PubMed PubMed Central

[22] Pang L., Shen L., Zhao Z., Mathematical modelling and analysis of the tumor treatment regimens with pulsed immunotherapy and chemotherapy, Comput Math Methods Med 2017, 6260474, 12.10.1155/2016/6260474Search in Google Scholar PubMed PubMed Central

[23] Polyanin A.D., Zaitsev V.F., Handbook of Exact Solutions for Ordinary Differential Equations 2003, Chapman and Hall/CRC.10.1201/9781420035339Search in Google Scholar

[24] Popovych R.O., More common errors in finding exact solutions of nonlinear differential equations: Part I, Commun Nonlinear Sci Numer Simul 2010, 15, 3887–3899.10.1016/j.cnsns.2010.01.037Search in Google Scholar

[25] Reluga T.C., Dahari H., Perelson A.S., Analysis of hepatitis c virus infection models with hepatocyte homeostasis, SIAM J Appl Math 2009, 69, 999–1023.10.1137/080714579Search in Google Scholar PubMed PubMed Central

[26] Scott A., Ed., Encyclopedia of Nonlinear Science 2004, Routledge, New York.Search in Google Scholar

[27] Wang J., Tian X., Global stability of a delay differential equation of hepatitis B virus infection with immune response, Electron J Differential Equations 2013, 1–11.Search in Google Scholar

[28] Wodarz D., Computational approaches to study oncolytic virus therapy: insights and challenges, Gene Ther Mol Biol 2004, 8, 137–146.Search in Google Scholar

[29] Wodarz D., Komarova N.L., Eds. Dynamics of Cancer 2014, World Scientific, Singapore.10.1142/8973Search in Google Scholar

Received: 2018-02-21

Accepted: 2018-11-23

Published Online: 2018-12-31

This work is licensed under the Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 License.

Articles in the same Issue

https://doi.org/10.1515/math-2018-0130

Keywords for this article

hepatitis C model; homoclinic/heteroclinic solution; generalized differential operator; inverse balancing

Creative Commons

BY-NC-ND 4.0