## 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.

## 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:

where p*t* is time; *β* 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; *d*T, *d*I are death rates for uninfected hepatocytes and infected cells respectively; *r*T, *r*I are parameters of the logistic proliferation of T and I respectively; logistic proliferation happens only if T < T* _{max}*; parameters

*η*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:

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:

where *c*, *u*, *v*, *a _{k}* ,

*b*P R,

_{k}*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*_{4} = *b*_{4} = 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:

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*_{0}, . . . , *a*_{4}; *b*_{0}, . . . , *b*_{4} 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:

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]:

where *η* ≠ 0, *σ*, γ ∈ R are constants; *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.

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

where *η* ≠ 0, *σ*, γ ∈ R are constants; *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:

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

where

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:

where *β*. Standard properties of differentiation operators hold true for (15) [14]:

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

#### 2.4.2 Multiplicative operators

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

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

Note that (20) follows immediately from the definition of (19). Without loss of generality, the proof of (21) for multiplicative operator

Let *y*_{1} :═ **M***u* ═ *y*_{1}(*t*, *u*, *v*), *y*_{2} :═ **M***v* ═ *y*_{2}(*t*, *u*, *v*), *z* :═ **M***f* (*u*, *v*) ═ *z*(*t*, *u*, *v*) and *w* :═ *f* (**M***u*,**M***v*) ═ *f* (*y*_{1}, *y*_{2}). To prove (21), it needs to be shown that *z* ═ *w* for all *t*, *u*, *v*.

Note that:

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

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

Analogously, it is shown that:

with *y*_{1}(0, *u*, *v*) ═ *u* and *y*_{2}(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*_{1}(0, *u*, *v*), *y*_{2}(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:

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

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

Substituting *t* for *t* ═ *c* yields (28).

#### 2.4.3 Construction of solutions to ODEs

Let us consider the following system of ODEs:

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

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

The convention **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

Analogously,

Selecting **G** yields that

In the following derivations,

will be used, which transforms (32) into:

Furthermore, coefficients *p _{j}* ,

*q*satisfy recurrence relations:

_{j}## 3 Existence of homoclinic/heteroclinic solutions in (30)

Let *ρ*_{1} ≠ *ρ*_{2}. If (30) admits solutions (13), (14) then (13) and (36) must be equal. Expanding (13) in a power series and equating to (36) yields:

Note that *p*_{0} ═ *u* by (35), thus (39) yields:

Analogous derivations with respect to *y* and *ν*_{1} ≠ *ν*_{2} result in:

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

**Theorem 3.1***System (30) admits homoclinic/heteroclinic solutions (13), (14) with ρ*_{1} ≠ *ρ*_{2}*if and only if*

*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:

Solving the above equations for *λ*_{1}, *λ*_{2} results in (44).

Equation (41) yields the following determinant equality:

Expanding the left side of (51) yields:

where

Solving (52) for *ρk* results in: _{k}

Since *ρ*_{1} ≠ *ρ*_{2}, the discriminant

Denoting

Using recursion (38) it can be obtained that:

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

Furthermore,

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

Applying operator

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

**Sufficiency** Condition (44) yields:

Applying operator

Continuing by induction yields (41).

The proof for parameters of *y* is analogous.

**Corollary 3.1***If conditions of Theorem 3.1 hold true, then the third and higher order Hankel determinants of sequences**are equal to zero*

*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:

with initial conditions

The following notations are introduced:

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

### 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:

Equation (78) results in:

Letting

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

thus (75) (and conversely (3)) only admits homoclinic/heteroclinic solutions with equal denominators. Let

where *T*(*t*) := (*t* = *(*_{1}) (*t* = *t*_{2}).

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

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

Note that

and

Taking *t* = *t*_{1}, *t*_{2} in (91) and using (93), (94) yields:

Analogous computations with respect to (92) result in:

Solution of (95)–(98) with respect to *a*_{2}, *a*_{3}, *b*_{2}, *b*_{3} reads:

Similarly, taking *a*_{0}, *a*_{4}, *b*_{0}, *b*_{4}:

Finally, taking *t* = 0 in (91), (92) yields *a*_{1}, *b*_{1}:

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*_{3} = *b*_{2} and *b*_{3} = *a*_{2} must hold if (75) admits solution (90).

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

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

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:

Parameter *η* must be chosen to satisfy

Furthermore, *η* can only depend on coefficients *a*_{0}, . . . , *a*_{4}; *b*_{0}, . . . , *b*_{4}, otherwise Theorem 3.1 does not hold true and obtained solutions would not be valid for all initial conditions.

It can be observed that:

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

and

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

then *A*_{6} = *A*_{4} = 0 and using (111), *η*^{2} can be expressed as:

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

where *α _{k}* ,

*β*are functions of

_{k}*a*

_{0}, . . . ,

*a*

_{4};

*b*

_{0}, . . . ,

*b*

_{4}.

Analogous computations with respect to

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

Note that:

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

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.1***System (3) admits homoclinic/heteroclinic solutions (8), (9) if and only if conditions (107), (108) and (123) hold true*.

Note that

where

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

where

Parameters *σ*,γ read:

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:

As demonstrated in [15], the value *η* 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:

To preserve biological significance of system (2), the parameters (134), (135) must satisfy *q*, *s*, *r* ≥ 0; *b* P [0^{-2}; 10^{3}]; *d* ∈ [10^{-3}; 10^{2}][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

Both equations are satisfied if parameter *r* reads:

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

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

**Corollary 6.1***Hepatitis C model (2) admits homoclinic/heteroclinic solutions if and only if**(137) and (138) hold true*.

Computer algebra computations prove that when Corollary (6.1) holds true, parameters *y*_{1} = *y*_{2} = 0.

### 6.2 Equilibria

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

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

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

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

### 6.3 Computational experiment

Let us consider the following system:

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

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:

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

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

where

Derivations given in Subsection 5.2 result in:

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

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:

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:

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

## 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*_{4}, *b*_{4} 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*_{4}, *b*_{4} → 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

© 2018 Telksnys et al., published by De Gruyter

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