Structural stability of the hepatitis C model with the proliferation of infected and uninfected hepatocytes

ABSTRACT

The structural stability of the hepatitis C model with the proliferation of infected and uninfected hepatocytes is investigated in this paper. The model is clinically verified to accurately reflect the viral dynamics of hepatitis C. It is demonstrated that the limit transition from the hepatitis C model to the Riccati system coupled with the diffusive terms leads to the elimination of the quadratic terms. Such a novel effect leads to the introduction of the concept of the deformed order-1 solitary solutions. The generalized operator of differentiation is used to construct the deformed order-1 solitary solutions to the Riccati system coupled with the diffusive terms. Finally, it is demonstrated that the Riccati system coupled with the diffusive terms admits non-deformed order-1 solitary solutions, which proves the structural instability of the hepatitis C model with the proliferation of infected and uninfected hepatocytes.

KEYWORDS:

• Hepatitis C model
• solitary solution
• structural stability

1. Introduction

The hepatitis C virus (HCV) infection represents a significant global public health problem [1]. One of the first mathematical models used to describe the kinetics of chronic HCV infection during treatment has been introduced more than two decades ago [2]. This paradigmatic phenomenological model include

s three ordinary coupled differential equations representing the population of target cells, productively-infected cells, and virus cells. However, the model in [2] is not able to explain some observed HCV kinetic profiles under treatment [3]. The model in [4] expands the HCV viral-dynamic model [2] by incorporating density-dependent proliferation what make the model predictions to agree with experimental observations during acute infection, under antiviral therapy, and after the cessation of therapy. The mathematical properties of this model, including steady state and dynamical behaviour are thoroughly analysed in [5].

A common feature of nonlinear dynamical systems is that solitary solutions represent a separatrix in the space of system parameters and initial conditions [6]. Therefore, the existence of solitary solutions may help understand the global dynamics of the system.

Kink solitary solutions to the hepatitis C model in [4] are derived in [7]. It is demonstrated in [7] that kink solutions can be in either linear or hyperbolic (inverse) relationship. If kink solutions are in the linear relationship, an infinitesimal perturbation of infected cell population results in an infinitesimal perturbation of uninfected cell population. However, if solutions are in the hyperbolic relationship, a small perturbation of the uninfected cell population can result in a large alteration of the infected cell population [7].

From a mathematical point of view, the hepatitis C model in [4] represents two nonlinear Riccati equations coupled with both multiplicative and diffusive terms [7]. Solitary solutions to Riccati system coupled with only multiplicative terms have already been studied in [8]. However, the dynamics of Riccati equations coupled with only diffusive terms remains unexplored territory.

Diffusively coupled models provide crucial insight into the dynamical behaviour of a wide range of biological and engineering systems. For example, synchronization of diffusively coupled models is a rich area of research [9]. In particular, the derivation of conditions that rule out synchronization does facilitate the study of spatial pattern formation [10,11].

Therefore, the main objective of this paper directly follows from the previous discussion. The existence of solitary solutions to a system of Riccati equations coupled with diffusive terms is an open question, particularly important due to the connections with the Hepatitis C model with the proliferation of infected and uninfected hepatocytes.

While diffusively coupled models are significant in modelling HCV evolution, numerous studies also consider their stochastic counterparts. A stochastic model of the spread of infectious diseases in considered in [12]. Numerous studies on stochastic models for the evolution of hepatitis B virus have been recently performed [13–16]. The impact of such models and some of their unifying properties is discussed in [17].

Another question of particular importance to the modelling of biological phenomena is boundedness of the solutions. This topic is discussed extensively in partial differential equation (PDE) models of biological systems, especially in parabolic PDEs. For example, unbounded solutions to a chemotaxis model are studied along with their blow-up criteria in [18] and [19].

The paper is structured as follows. Basic definitions are given in the Preliminaries section. The concept of the deformed order-1 solitary solutions is presented in section 3. The role of the deformation function is discussed in section 4. The structure and the simplification of the deformed system is derived in sections 5 and 6. Deformed order-1 solitary solutions are constructed in section 6. The structural stability of the hepatitis C model with the proliferation of infected and uninfected hepatocytes is discussed in section 7. Concluding remarks are given in the last section.

2. Preliminaries

2.1. The hepatitis C virus infection model

The hepatitis C virus infection model with the inclusion of the proliferation of infected and uninfected hepatocytes reads [5]:

(1) $\left\{\begin{array}{c}\frac{\mathit{dT}}{\mathit{dt}}=s_0+r_1\mathit{T}\left(1-\frac{\mathit{T}+\mathit{I}}{{\mathit{T}}_{\mathit{max}}}\right)-d_1\mathit{T}-(1-\mathit{\eta })\mathit{\beta VT}+\mathit{qI},\\ \frac{\mathit{dI}}{\mathit{dt}}=r_2\mathit{I}\left(1-\frac{\mathit{T}+\mathit{I}}{{\mathit{T}}_{\mathit{max}}}\right)-d_2\mathit{I}+(1-\mathit{\eta })\mathit{\beta VT}-\mathit{qI},\\ \frac{\mathit{dV}}{\mathit{dt}}=(1-\mathit{\epsilon })\mathit{pI}-\mathit{cV},\end{array}\right.$(1)

where $\mathit{t}$ is time, $\mathit{T}(\mathit{t})$ represents uninfected hepatocytes; $\mathit{I}(\mathit{t})$ represents infected cells; $\mathit{V}(\mathit{t})$ represents virus-free population. All other model parameters are constants: $\mathit{s}$ and $\mathit{q}$ represent the increase rate of uninfected hepatocytes through immigration and spontaneous cure by noncytolytic process respectively; $r_1$ and $r_2$ are parameters of the logistic proliferation of $\mathit{T}$ and $\mathit{I}$ respectively (logistic proliferation happens only if $\mathit{T}<T_{\mathit{max}}$); $d_1$ and $d_2$ are death rates for uninfected hepatocytes and infected cells respectively; $\mathit{\eta }$ is the infection reduction rate induced by the antiviral treatment; $\mathit{\beta }$ is the rate of infection per free virus per hepatocyte; $\mathit{\epsilon }$ is the fraction of the viral production rate induced by the antiviral treatment; $\mathit{p}$ is the free virus production rate per infected cell; $\mathit{c}$ is the immune virus clearance rate.

It is demonstrated in [5] that the range of rates of viral clearance is significantly faster compared to other parameters on the time scale. Furthermore, the viral dynamics closely follow the dynamics of infected after a short transient process, resulting in the following equality for $\mathit{V}(\mathit{t})$ [5]:

(2) $\mathit{V}(\mathit{t})=\left(1-\mathit{\epsilon }\right)\mathit{V}\left(t_0\right)+\mathit{\epsilon V}\left(t_0\right)exp\left(-\mathit{c}\left(\mathit{t}-t_0\right)\right).$(2)

Then, (1) can be rearranged into the general form by the introduction of dimensionless state variables [7]:

(3) $\left\{\begin{array}{c}\frac{\mathit{dx}}{\mathit{dt}}=a_0+a_1\mathit{x}+a_2{x}^2+a_3\mathit{xy}+a_4\mathit{y},\\ \frac{\mathit{dy}}{\mathit{dt}}=b_0+b_1\mathit{y}+b_2{y}^2+b_3\mathit{xy}+b_4\mathit{x}.\end{array}\right.$(3)

where $a_k,b_k\in \mathbb{R}$; $\mathit{k}=0,1,\dots ,4$. From the mathematical point of view, system (3) does represent two Riccati differential equations with constant coefficients coupled with diffusive and multiplicative terms.

It is worthwhile to note that (3) represents a classical system for modelling HCV dynamics. These systems may have a variety of modifications, such as the application of a fractional-order derivative instead of classical derivation. The properties of such of such derivatives have been discussed in detail in [20]. The value of their application to biological systems has been established in [21] and a particular case concerning the hepatitis B virus (HBV) is considered in [22].

2.2. Solitary solutions and their orders

The standard form of a solitary solution reads [23]:

(4) $\mathit{x}\left(\mathit{t}\right)=\mathit{\sigma }\frac{{\prod }_{\mathit{k}=1}^n\left(exp\left(\mathit{\eta }\left(\mathit{t}-t_0\right)\right)-x_k\right)}{{\prod }_{\mathit{k}=1}^n\left(exp\left(\mathit{\eta }\left(\mathit{t}-t_0\right)\right)-t_k\right)},$(4)

where $\mathit{n}\in \mathbb{N}$ is the order of the solitary solution; and $\mathit{\sigma }$, $\mathit{\eta }$, $t_0$, $x_k$, $t_k\in \mathbb{R}$ are real parameters.

The first order solitary solution (at $\mathit{n}=1$) represents a sigmoid function describing the transition from one steady state to another steady state via a monotonous trajectory.

2.3. Riccati equation with constant coefficients

The Riccati equation with constant coefficients

(5) $\frac{\mathit{dx}}{\mathit{dt}}=a_0+a_1\mathit{x}+a_2{x}^2$(5)

admits the first order solitary solution when parameters $a_0$, $a_1$, $a_2\in \mathbb{R}$ are such that the roots of the polynomial $a_0+a_1\mathit{x}+a_2{x}^2=a_2\left(\mathit{x}-r_1\right)\left(\mathit{x}-r_2\right)$ are real $\left(r_1\le r_2\right)$. Then, the first order solitary solution to 5 reads [6]:

(6) $\mathit{x}(\mathit{t})=\mathit{\sigma }\frac{exp\left(\eta \left(t-t_0\right)\right)-x_1}{exp\left(\mathit{\eta }\left(\mathit{t}-t_0\right)\right)-t_1}$(6)

where $\mathit{\sigma }=r_2$; $\mathit{\eta }=a_2\left(r_1-r_2\right)$; $x_1=\frac{r_1}{r_2}\cdot \frac{u-r_2}{\mathit{u}-r_1}$; $t_1=\frac{u-r_2}{\mathit{u}-r_1}$; $\mathit{u}$ denotes the initial condition $\mathit{x}\left(t_0\right)=\mathit{u}$ satisfying the system of inequalities $r_1\le \mathit{u}\le r_2$. The first order solitary solution is depicted in Figure 1.

Figure 1. The first order solitary solution to (5) at $a_2=1,r_1=\frac{1}{2},r_2=-\frac{2}{3}$ and the initial condition $\mathit{x}(0)=\frac{1}{3}$.

2.4. A system of diffusively coupled Riccati equations

The paradigmatic model of two diffusively coupled dynamical systems reads [9]:

(7) $\left\{\begin{array}{c}\frac{\mathit{dx}}{\mathit{dt}}=\mathit{F}\left(\mathit{x}\right)+{\epsilon }_x\left(\mathit{y}-\mathit{x}\right);\\ \frac{\mathit{dy}}{\mathit{dt}}=\mathit{G}\left(\mathit{y}\right)+{\epsilon }_y\left(\mathit{x}-\mathit{y}\right),\end{array}\right.$(7)

where $\mathit{t}$ is time; $\mathit{F}$ and $\mathit{G}$ are the functions that represent isolated chaotic dynamics of $\mathit{x}$ and $\mathit{y}$; ${\epsilon }_x$ and ${\epsilon }_y$ are diffusive coupling parameters usually set as a positive constants.

Thus, the system of two diffusively coupled Riccati equations can be reduced to the following standard form:

(8) $\left\{\begin{array}{c}\frac{\mathit{dx}}{\mathit{dt}}=a_0+a_1\mathit{x}+a_2{x}^2+{\epsilon }_x\mathit{y};\\ \frac{\mathit{dy}}{\mathit{dt}}=b_0+b_1\mathit{y}+b_2{y}^2+{\epsilon }_y\mathit{x},\end{array}\right.$(8)

where $a_k,b_k\in \mathbb{R}$, $\mathit{k}=1,2,3$; $a_2\ne 0$; $b_2\ne 0$.

3. The concept of the deformed order–1 solitary solutions to (8)

Definition 3.1.

The order–1 solitary solutions to 8 read:

(9) $\left\{\begin{array}{c}\mathit{x}(\mathit{t})={\sigma }_x\frac{exp\left(\eta \left(t-t_0\right)\right)-x_1}{exp\left(\mathit{\eta }\left(\mathit{t}-t_0\right)\right)-t_1};\\ \mathit{y}(\mathit{t})={\sigma }_y\frac{exp\left(\eta \left(t-t_0\right)\right)-y_1}{exp\left(\mathit{\eta }\left(\mathit{t}-t_0\right)\right)-t_1}.\end{array}\right.$(9)

Note that parameters $\mathit{\eta }$, $t_0$, and $t_1$ must be the same for both functions $\mathit{x}(\mathit{t})$ and $\mathit{y}(\mathit{t})$ in order to ensure the preservation of consistency between the coupled nonlinear differential equations [7].

3.1. The existence of the order–1 solitary solutions to (8)

The existence of the order–1 solitary solutions to (8) follows directly from the properties of (3). The necessary condition for the existence of the order–1 solitary solutions to (3) is given by [7]:

(10) $\begin{array}{cc}{a}_3& =b_2;\\ b_3& =a_2.\end{array}$(10)

The model (3) is mapped to (8) when both parameters $a_3$ and $b_3$ are set to zero (the multiplicative coupling is eliminated from (3)). However, in that case, the conditions for the existence of the order–1 solitary solutions (10) require that $a_2=0$ and $b_2=0$ which contradicts the definition of the Riccati equations (5)(5) $\frac{\mathit{dx}}{\mathit{dt}}=a_0+a_1\mathit{x}+a_2{x}^2$(5) (Figure 2). Therefore, the following hypothesis is posed.

Figure 2. A schematic diagram illustrating the limit transition from the hepatitis C model to the degenerate riccati system with diffusive coupling.

Hypothesis 3.1.

The limit transition from (3) to (8) suggests that order–1 solitary solutions to (8) do not exist.

3.2. The proposed architecture of the deformed order–1 solitary solutions

Definition 3.2.The deformed order–1 solitary solution reads:

(11) $\left\{\begin{array}{c}\mathit{x}(\mathit{t})={\sigma }_x\frac{\phi (t)-x_1}{\mathit{\phi }(\mathit{t})-t_1};\\ \mathit{y}(\mathit{t})={\sigma }_y\frac{\phi (t)-y_1}{\mathit{\phi }(\mathit{t})-t_1},\end{array}\right.$(11)

where the function $\mathit{\phi }(\mathit{t})$ defines the unknown time deformation.

In other words, a hypothesis is posed that (8) admits the order–1 solitary solutions if only the time axis is deformed according to the deformation function $\mathit{\phi }(\mathit{t})$.

Remark 1. For any pair of deformed order-1 solitary solutions (11), there exists parameters $\mathit{A},\mathit{B}\in \mathbb{R}$ such that $\mathit{x}(\mathit{t})$ and $\mathit{y}(\mathit{t})$ are in a linear relationship $\mathit{Ax}(\mathit{t})+\mathit{By}(\mathit{t})=1$.

Proof. Using basic rearrangements, (11) yields:

(12) ${\sigma }_y\left(t_1-y_1\right)\mathit{x}(\mathit{t})-{\sigma }_x\left(t_1-x_1\right)\mathit{y}(\mathit{t})={\sigma }_x{\sigma }_y\left(x_1-y_1\right).$(12)

Dividing both sides by ${\sigma }_x{\sigma }_y\left(x_1-y_1\right)$ yields:

(13) $\mathit{A}=\frac{t_1-y_1}{{\sigma }_x\left(x_1-y_1\right)},\mathit{B}=\frac{t_1-x_1}{{\sigma }_y\left(y_1-x_1\right)},$(13)

which proves the statement.

4. The role of the deformation function $\mathit{\phi }(\mathit{t})$

The proposed deformed order–1 solitary solution (11) may not satisfy (8). Note that (11) defines the non-deformed order–1 solitary solution (9) at $\mathit{\phi }(\mathit{t})=exp\mathit{\eta }(\mathit{t}-t_0)$. Therefore, (8) is extended in order to accommodate a larger class of solutions:

(14) $\left\{\begin{array}{c}\mathit{f}(\mathit{t})\frac{\mathit{dx}}{\mathit{dt}}=a_0+a_1\mathit{x}+a_2{x}^2+{\epsilon }_x\mathit{y},\\ \mathit{f}(\mathit{t})\frac{\mathit{dy}}{\mathit{dt}}=b_0+b_1\mathit{y}+b_2{y}^2+{\epsilon }_y\mathit{x},\end{array}\right.$(14)

where $\mathit{f}(\mathit{t})$ is an unknown function representing this extension.

Definition 4.1. System (14) is called the extended system of (8).

Definition 4.2. The transformed time scale is defined as $\hat{t}:=\mathit{\phi }(\mathit{t})$.

The variable substitution ($\mathit{t}$ by $\hat{t}$) modifies the structure of (8). These changes are represented by the function $\hat{f}\left(\hat{t}\right)$:

(15) $\left\{\begin{array}{c}\hat{f}\left(\hat{t}\right)\frac{d\hat{x}}{\mathit{d}\hat{t}}=a_0+a_1\hat{x}+a_2{\hat{x}}^2+{\epsilon }_x\hat{y},\\ \\ \hat{f}\left(\hat{t}\right)\frac{d\hat{y}}{\mathit{d}\hat{t}}=b_0+b_1\hat{y}+b_2{\hat{y}}^2+{\epsilon }_y\hat{x}.\end{array}\right.$(15)

Definition 4.3. System (15) is called the image system of (14).

Lemma 4.4. Functions $\hat{f}\left(\hat{t}\right)$, $\hat{t}=\mathit{\phi }(\mathit{t})$, and $\mathit{f}(\mathit{t})$ (if they exist) are related by the following differential equation:

(16) $\frac{\mathit{d\phi }(\mathit{t})}{\mathit{dt}}=\frac{\hat{f}(\mathit{\phi }(\mathit{t}))}{\mathit{f}(\mathit{t})}.$(16)

Proof. Consider the following solution to the extended system (14): $\mathit{x}=\mathit{x}(\mathit{t})$; $\mathit{y}=\mathit{y}(\mathit{t})$. Let us assume that the image system (15) exists. Then the solution to the image system (15) does satisfy the following identities:

(17) $\begin{array}{cc}\mathit{x}(\mathit{t})& =\hat{x}(\mathit{\phi }(\mathit{t}));\\ \mathit{y}(\mathit{t})& =\hat{y}(\mathit{\phi }(\mathit{t})),\end{array}$(17)

if only the variable change function $\hat{t}:=\mathit{\phi }(\mathit{t})$ exists. Then,

(18) $\left\{\begin{array}{c}\frac{\mathit{dx}(\mathit{t})}{\mathit{dt}}=\frac{d\hat{x}\left(\hat{t}\right)}{\mathit{d}\hat{t}}{|}_{\hat{t}=\mathit{\phi }(\mathit{t})}\cdot \frac{\mathit{d\phi }(\mathit{t})}{\mathit{dt}};\\ \frac{\mathit{dy}(\mathit{t})}{\mathit{dt}}=\frac{d\hat{y}\left(\hat{t}\right)}{\mathit{d}\hat{t}}{|}_{\hat{t}=\mathit{\phi }(\mathit{t})}\cdot \frac{\mathit{d\phi }(\mathit{t})}{\mathit{dt}}.\end{array}\right.$(18)

The image system can be transformed into:

(19) $\left\{\begin{array}{c}\left(\hat{f}\left(\hat{t}\right)\frac{d\hat{x}}{\mathit{d}\hat{t}}\right){|}_{\hat{t}=\mathit{\phi }(\mathit{t})}=\left(a_0+a_1\hat{x}+a_2{\hat{x}}^2+{\epsilon }_x\hat{y}\right){|}_{\hat{t}=\mathit{\phi }(\mathit{t})},\\ \left(\hat{f}\left(\hat{t}\right)\frac{d\hat{y}}{\mathit{d}\hat{t}}\right){|}_{\hat{t}=\mathit{\phi }(\mathit{t})}=\left(b_0+b_1\hat{y}+b_2{\hat{y}}^2+{\epsilon }_y\hat{x}\right){|}_{\hat{t}=\mathit{\phi }(\mathit{t})}.\end{array}\right.$(19)

Therefore,

(20) $\left\{\begin{array}{c}\hat{f}(\mathit{\phi }(\mathit{t}))\cdot \frac{\mathit{dx}(\mathit{t})}{\mathit{dt}}/\frac{\mathit{d\phi }(\mathit{t})}{\mathit{dt}}=a_0+a_1\mathit{x}+a_2{x}^2+{\epsilon }_x\mathit{y},\\ \hat{f}(\mathit{\phi }(\mathit{t}))\cdot \frac{\mathit{dy}(\mathit{t})}{\mathit{dt}}/\frac{\mathit{d\phi }(\mathit{t})}{\mathit{dt}}=b_0+b_1\mathit{y}+b_2{y}^2+{\epsilon }_y\mathit{x}.\end{array}\right.$(20)

Comparing (14) and (20) yields the following identities:

(21) $\left\{\begin{array}{c}\hat{f}(\mathit{\phi }(\mathit{t}))\cdot \frac{\mathit{dx}(\mathit{t})}{\mathit{dt}}/\frac{\mathit{d\phi }(\mathit{t})}{\mathit{dt}}=\mathit{f}(\mathit{t})\frac{\mathit{dx}(\mathit{t})}{\mathit{dt}};\\ \hat{f}(\mathit{\phi }(\mathit{t}))\cdot \frac{\mathit{dy}(\mathit{t})}{\mathit{dt}}/\frac{\mathit{d\phi }(\mathit{t})}{\mathit{dt}}=\mathit{f}(\mathit{t})\frac{\mathit{dy}(\mathit{t})}{\mathit{dt}},\end{array}\right.$(21)

which concludes the proof.

Corollary 4.5 $\mathit{x}(\mathit{t})=\hat{x}\left(\mathit{\phi }(\mathit{t})\right)$; $\mathit{y}(\mathit{t})=\hat{y}\left(\mathit{\phi }(\mathit{t})\right)$.

Example 4.6. The Exp-function method for the Riccati equation.

The Exp-function method, which was proposed more than two decades ago in [24], uses the exp-function variable substitution in order to transform the original differential equation to the image differential equation. Let us consider the Riccati equation with constant coefficients (5). The new variable $\hat{t}$ is set in the form of the exponent:

(22) $\hat{t}=exp\mathit{\eta }(\mathit{t}-t_0)=\mathit{\phi }(\mathit{t}).$(22)

The change of variables transforms (5) into the image differential equation:

(23) $\mathit{\eta }\hat{t}\frac{d\hat{x}}{\mathit{d}\hat{t}}=a_0+a_1\hat{x}+a_2{\hat{x}}^2.$(23)

The solution to (23) can be constructed using operator techniques [25] and expressed in the closed form:

(24) $\mathit{x}\left(\hat{t}\right)=r_2\frac{\hat{t}-x_1}{\hat{t}-t_1},$(24)

where $r_2$, $x_1$, and $t_1$ are defined in (6). Therefore, $\hat{f}\left(\hat{t}\right)=\mathit{\eta }\hat{t}$ and $\mathit{f}(\mathit{t})=1$. It follows from (22) that $\frac{\mathit{d\phi }}{\mathit{dt}}=\mathit{\eta }exp\mathit{\eta }(\mathit{t}-t_0)$. On the other hand, according to (16):

(25) $\frac{\mathit{d\phi }}{\mathit{dt}}=\frac{\hat{f}\left(\hat{t}\right)}{\mathit{f}(\mathit{t})}=\mathit{\eta }\hat{t}=\mathit{\eta }exp\mathit{\eta }(\mathit{t}-t_0).$(25)

Thus, the introduction of the deformation function $\mathit{\phi }(\mathit{t})$ can be considered as the generalization of the Exp-function method where the function $\mathit{f}(\mathit{t})$ is not necessarily set to $1$.

Example 4.7. The classical Exp-function method implies $\hat{t}=exp(\mathit{\eta t})$. Then, functions $\hat{f}\left(\hat{t}\right)=\mathit{\eta }\hat{t}$, $\hat{t}=exp(\mathit{\eta t})$ and $\mathit{f}(\mathit{t})=1$ are related by (16). Thus, the image system (15) is given by:

(26) $\left\{\begin{array}{c}\mathit{\eta }\hat{t}\frac{d\hat{x}}{\mathit{d}\hat{t}}=a_0+a_1\hat{x}+a_2{\hat{x}}^2+{\epsilon }_x\hat{y},\\ \\ \mathit{\eta }\hat{t}\frac{d\hat{y}}{\mathit{d}\hat{t}}=b_0+b_1\hat{y}+b_2{\hat{y}}^2+{\epsilon }_y\hat{x},\end{array}\right.$(26)

and $\hat{x}={\sigma }_x\frac{\hat{t}-x_1}{\hat{t}-t_1}$; $\hat{y}={\sigma }_y\frac{\hat{t}-y_1}{\hat{t}-t_1}$ are the solutions to (26).

In other words, the original problem is split into two problems. The process of the construction of the solution to (14) is transformed into two consecutive problems (16) and (15).

5. The structure of the function $\hat{f}\left(\hat{t}\right)$

As previously demonstrated, the deformed order–1 solitary solution to the image system (15) reads:

(27) $\left\{\begin{array}{c}\hat{x}\left(\hat{t}\right)={\sigma }_x\frac{\hat{t}-x_1}{\hat{t}-t_1},\\ \hat{y}\left(\hat{t}\right)={\sigma }_y\frac{\hat{t}-y_1}{\hat{t}-t_1}.\end{array}\right.$(27)

Thus,

(28) $\frac{d\hat{x}\left(\hat{t}\right)}{\mathit{d}\hat{t}}={\sigma }_x\frac{x_1-t_1}{{(\hat{t}-t_1)}^2},\frac{d\hat{y}\left(\hat{t}\right)}{\mathit{d}\hat{t}}={\sigma }_y\frac{y_1-t_1}{{(\hat{t}-t_1)}^2}.$(28)

Inserting (27) and (28) into (15) yields a system of two identities:

(29) $\left\{\begin{array}{c}\begin{array}{cc}\hat{f}\left(\hat{t}\right){\sigma }_x(x_1-t_1)& =a_0{(\hat{t}-t_1)}^2+a_1{\sigma }_x(\hat{t}-t_1)(\hat{t}-x_1)+a_2{\sigma }_x^2{(\hat{t}-x_1)}^2\\ \mathit{\hspace{0.17em}}& +{\epsilon }_x{\sigma }_y(\hat{t}-y_1)(\hat{t}-t_1),\end{array}\\ \begin{array}{cc}\hat{f}\left(\hat{t}\right){\sigma }_y(y_1-t_1)& =b_0{(\hat{t}-t_1)}^2+b_1{\sigma }_y(\hat{t}-t_1)(\hat{t}-y_1)+b_2{\sigma }_y^2{(\hat{t}-y_1)}^2\\ \mathit{\hspace{0.17em}}& +{\epsilon }_y{\sigma }_x(\hat{t}-x_1)(\hat{t}-t_1).\end{array}\end{array}\right.$(29)

Definition 5.1. System (29) is called the balancing system of (15).

Taking $\hat{t}$ equal to $t_1$, $x_1$, and $y_1$ in the balancing system (29) yields the following systems of equations:

(30) $\left\{\begin{array}{c}\hat{f}(t_1){\sigma }_x(x_1-t_1)=a_2{\sigma }_x^2{(t_1-x_1)}^2,\\ \hat{f}(t_1){\sigma }_y(y_1-t_1)=b_2{\sigma }_y^2{(t_1-y_1)}^2;\end{array}\right.$(30)

(31) $\left\{\begin{array}{c}\hat{f}(x_1){\sigma }_x(x_1-t_1)=a_0{(x_1-t_1)}^2+{\epsilon }_x{\sigma }_y(x_1-y_1)(x_1-t_1),\\ \begin{array}{cc}\hat{f}(x_1){\sigma }_y(y_1-t_1)& =b_0{(x_1-t_1)}^2+b_1{\sigma }_y(x_1-y_1)(x_1-t_1)\\ \mathit{\hspace{0.17em}}& +b_2{\sigma }_y^2{(x_1-y_1)}^2;\end{array}\end{array}\right.$(31)

(32) $\left\{\begin{array}{c}\begin{array}{cc}\hat{f}(y_1){\sigma }_x(x_1-t_1)& =a_0{(y_1-t_1)}^2+a_1{\sigma }_x(y_1-x_1)(y_1-t_1)\\ \mathit{\hspace{0.17em}}& +a_2{\sigma }_x^2{(y_1-x_1)}^2,\end{array}\\ \hat{f}(y_1){\sigma }_y(y_1-t_1)=b_0{(y_1-t_1)}^2+{\epsilon }_y{\sigma }_x(y_1-x_1)(y_1-t_1);\end{array}\right.$(32)

Note that $\hat{f}(t_1)$, $\hat{f}(x_1)$, and $\hat{f}(y_1)$ are constants (the values of the function $\hat{f}\left(\hat{t}\right)$ at $\hat{t}=t_1$, $x_1$, and $y_1$ accordingly).

In total, Riccati equations in (15) are comprised of eight unknowns ($a_k,b_k,\mathit{k}=0,1,2$ and ${\epsilon }_x,{\epsilon }_y$). Computer algebra helps to solve the system for six unknowns $a_0$, $a_1$, $a_2$, $b_0$, $b_1$, and $b_2$ from ((30), (31), (32)):

(33) $a_0=-\frac{\hat{f}(x_1){\sigma }_x-(x_1-y_1){\epsilon }_x{\sigma }_y}{t_1-x_1}$(33)

(34) $\begin{array}{cc}{a}_1& =\frac{\hat{f}(x_1)(t_1-y_1{)}^2{\sigma }_x-\hat{f}(y_1)(t_1-x_1{)}^2{\sigma }_x+\hat{f}(t_1)(x_1-y_1{)}^2{\sigma }_x}{(t_1-x_1)(t_1-y_1)(x_1-y_1){\sigma }_x}\\ \mathit{\hspace{0.17em}}& -\frac{{(t_1-y_1)}^2(x_1-y_1){\epsilon }_x{\sigma }_y}{(t_1-x_1)(t_1-y_1)(x_1-y_1){\sigma }_x}\end{array}$(34)

(35) $a_2=-\frac{\hat{f}(t_1)}{(t_1-x_1){\sigma }_x}$(35)

(36) $b_0=-\frac{\hat{f}(y_1){\sigma }_y-(y_1-x_1){\epsilon }_y{\sigma }_x}{t_1-y_1}$(36)

(37) $\begin{array}{cc}{b}_1& =\frac{\hat{f}(y_1)(t_1-x_1{)}^2{\sigma }_y-\hat{f}(x_1)(t_1-y_1{)}^2{\sigma }_y+\hat{f}(t_1)(y_1-x_1{)}^2{\sigma }_y}{(t_1-x_1)(t_1-y_1)(y_1-x_1){\sigma }_y}\\ \mathit{\hspace{0.17em}}& -\frac{{(t_1-x_1)}^2(y_1-x_1){\epsilon }_y{\sigma }_x}{(t_1-x_1)(t_1-y_1)(y_1-x_1){\sigma }_y}\end{array}$(37)

(38) $b_2=-\frac{\hat{f}(t_1)}{(t_1-y_1){\sigma }_y}$(38)

Inserting the expressions of $a_0$, $a_1$, $a_2$, $b_0$, $b_1$, and $b_2$ into any equality of the balancing system (29) yields the same expression of $\hat{f}\left(\hat{t}\right)$:

(39) $\begin{array}{cc}\hat{f}\left(\hat{t}\right)& =\frac{-\hat{f}(x_1)(t_1-y_1)(\hat{t}-t_1)(\hat{t}-y_1)+\hat{f}(y_1)(t_1-x_1)(\hat{t}-t_1)(\hat{t}-x_1)}{(t_1-x_1)(t_1-y_1)(x_1-y_1)}\\ \mathit{\hspace{0.17em}}& +\frac{\hat{f}(t_1)(x_1-y_1)(\hat{t}-x_1)(\hat{t}-y_1)}{(t_1-x_1)(t_1-y_1)(x_1-y_1)}.\end{array}$(39)

It is interesting to note that $\hat{f}\left(\hat{t}\right)$ is not dependant on ${\epsilon }_x,{\epsilon }_y$. Moreover, $\hat{f}\left(\hat{t}\right)$ is a second order polynomial with respect of $\hat{t}$:

(40) $\hat{f}\left(\hat{t}\right)=A_0+A_1\hat{t}+A_2{\hat{t}}^2,$(40)

where

(41) $A_0=\frac{-\hat{f}(x_1)t_1{y}_1(t_1-y_1)+\hat{f}(y_1)t_1{x}_1(t_1-x_1)+\hat{f}(t_1)x_1{y}_1(x_1-y_1)}{(t_1-x_1)(t_1-y_1)(x_1-y_1)},$(41)

(42) $A_1=\frac{\hat{f}(x_1)(t_1^2-y_1^2)-\hat{f}(y_1)(t_1^2-x_1^2)-\hat{f}(t_1)(x_1^2-y_1^2)}{(t_1-x_1)(t_1-y_1)(x_1-y_1)},$(42)

(43) $A_2=\frac{-\hat{f}(x_1)(t_1-y_1)+\hat{f}(y_1)(t_1-x_1)+\hat{f}(t_1)(x_1-y_1)}{(t_1-x_1)(t_1-y_1)(x_1-y_1)}.$(43)

6. The simplification of the structure of the function $\hat{f}\left(\hat{t}\right)$

The balancing system (29) reveals that the coefficient $A_2$ in (40) is equal to zero if the following identities hold true:

(44) $\left\{\begin{array}{c}{a}_0+a_1{\sigma }_x+a_2{\sigma }_x^2+{\epsilon }_x{\sigma }_y=0,\\ b_0+b_1{\sigma }_y+b_2{\sigma }_y^2+{\epsilon }_y{\sigma }_x=0.\end{array}\right.$(44)

Let us assume that the parameters $a_0$, $a_1$, $a_2$, $b_0$, $b_1$, $b_2$, ${\epsilon }_x$, and ${\epsilon }_y$ are given by the model of the system of two Riccati equations coupled with diffusive terms (8). Then, the parameters ${\sigma }_x$ and ${\sigma }_y$ can be computed via (44). In other words, a proper selection of ${\sigma }_x$ and ${\sigma }_y$ (the parameters of the deformed order–1 solitary solution) ensures that identities (44) hold true.

Lemma 6.1. Let the parameters $a_0$, $a_1$, $a_2$, $b_0$, $b_1$, $b_2$, ${\epsilon }_x$, ${\epsilon }_y$, $t_1$ be fixed, and the parameters ${\sigma }_x$, ${\sigma }_y$, $x_1$, $x_2$ satisfy (44) in addition to the following system of identities:

(45) $\left\{\begin{array}{c}{\sigma }_y(y_1-t_1)(a_0{t}_1^2+a_1{\sigma }_x{t}_1{x}_1+a_2{\sigma }_x^2{x}_1^2+{\epsilon }_x{\sigma }_y{y}_1{t}_1)=\\ ={\sigma }_x(x_1-t_1)(b_0{t}_1^2+b_1{\sigma }_y{t}_1{y}_1+b_2{\sigma }_y^2{y}_1^2+{\epsilon }_y{\sigma }_x{x}_1{t}_1),\\ {\sigma }_y(y_1-t_1)(2a_0{t}_1+a_1{\sigma }_x(t_1+x_1)+2a_2{\sigma }_x^2{x}_1+{\epsilon }_x{\sigma }_y(y_1+t_1))=\\ ={\sigma }_x(x_1-t_1)(2b_0{t}_1+b_1{\sigma }_y(t_1+y_1)+2b_2{\sigma }_y^2{y}_1+{\epsilon }_y{\sigma }_x(x_1+t_1)).\end{array}\right.$(45)

Then, the structure of the function $\hat{f}\left(\hat{t}\right)$ is simplified to:

(46) $\hat{f}\left(\hat{t}\right)=A_0+A_1\hat{t},$(46)

and the deformed order–1 solitary solution (27) satisfies the image system (15).

Proof. The proof follows directly from Lemma 4.4 and the balancing system (29).

Taking Lemma 6.1 into account, expressions (33)–(38) are simplified into:

(47) $a_0=\frac{A_0{\sigma }_x+A_1{\sigma }_x{x}_1-{\sigma }_y{x}_1{\epsilon }_x+{\sigma }_y{y}_1{\epsilon }_x}{x_1-t_1};$(47)

(48) $a_1=\frac{2A_0{\sigma }_x+A_1{\sigma }_x{t}_1+A_1{\sigma }_x{x}_1-{\sigma }_y{t}_1{\epsilon }_x+{\sigma }_y{y}_1{\epsilon }_x}{{\sigma }_x\left(t_1-x_1\right)};$(48)

(49) $a_2=\frac{A_0+A_1{t}_1}{{\sigma }_x\left(x_1-t_1\right)};$(49)

(50) $b_0=\frac{A_0{\sigma }_y+A_1{\sigma }_y{y}_1-{\sigma }_x{y}_1{\epsilon }_y+{\sigma }_x{x}_1{\epsilon }_y}{y_1-t_1};$(50)

(51) $b_1=\frac{2A_0{\sigma }_y+A_1{\sigma }_y{t}_1+A_1{\sigma }_y{y}_1-{\sigma }_x{t}_1{\epsilon }_y+{\sigma }_x{x}_1{\epsilon }_y}{{\sigma }_y\left(t_1-y_1\right)};$(51)

(52) $b_2=\frac{A_0+A_1{t}_1}{{\sigma }_y\left(y_1-t_1\right)}.$(52)

Example 6.2. This example illustrates the inverse balancing technique used to reconstruct the extended system of two Riccati equations coupled with diffusive terms (14) from the deformed order–1 solitary solution (27).

Consider parameter values ${\sigma }_x=-\frac{7}{4}$, ${\sigma }_y=\frac{19}{6}$, $t_1=-2$, $x_1=-\frac{22}{7}$, $y_1=-\frac{50}{19}$. Thus,

(53) $\left\{\begin{array}{c}\hat{x}\left(\hat{t}\right)=-\frac{22+7\hat{t}}{8+4\hat{t}},\\ \hat{y}\left(\hat{t}\right)=\frac{50+19\hat{t}}{12+6\hat{t}}.\end{array}\right.$(53)

Let us choose $\hat{f}(t_1)=-2$, $\hat{f}\left(x_1\right)=-\frac{22}{7}$, and $\hat{f}(y_1)=-\frac{50}{19}$ which guarantees the identity $-\hat{f}(x_1)(t_1-y_1)+\hat{f}(y_1)(t_1-x_1)+\hat{f}(x_1-y_1)=0$. Therefore, by (43) $A_2=0$ ($A_0=0$, $A_1=1$).

Now, ${\epsilon }_x=\frac{1}{2}$ and ${\epsilon }_y=\frac{1}{3}$ yield coefficients $a_0=-\frac{265}{48}$, $a_1=-4$, $a_2=-1$, $b_0=\frac{229}{18}$, $b_1=-7$, $b_2=1$.

Let us investigate two cases of function $\mathit{\phi }(\mathit{t})$:

Case 1. Let $\mathit{\phi }(\mathit{t})=t^2$. Then, $\hat{f}\left(\hat{t}\right)=\hat{t}$ and $\mathit{f}(\mathit{t})=\frac{t}{2}$. Finally,

(54) $\left\{\begin{array}{c}\mathit{x}(\mathit{t})=-\frac{22+7t^2}{8+4t^2},\\ \mathit{y}(\mathit{t})=\frac{50+19t^2}{12+6t^2}\end{array}\right.$(54)

is the analytical solution (see Figure 3) to:

(55) $\left\{\begin{array}{c}\frac{t}{2}\frac{\mathit{dx}(\mathit{t})}{\mathit{dt}}=-\frac{265}{48}-4\mathit{x}-x^2+\frac{1}{2}\mathit{y},\\ \frac{t}{2}\frac{\mathit{dy}(\mathit{t})}{\mathit{dt}}=\frac{229}{18}-7\mathit{y}+y^2+\frac{1}{3}\mathit{x}.\end{array}\right.$(55)

Figure 3. Deformed order-1 solitary solutions $\mathit{x}(\mathit{t}),\mathit{y}(\mathit{t})$ (see eq. (54)) to (55).

Case 2. Let $\mathit{\phi }(\mathit{t})=sin(\mathit{t})$. Then, $\hat{f}\left(\hat{t}\right)=\hat{t}$ and $\mathit{f}(\mathit{t})=tan(\mathit{t})$. Finally,

(56) $\left\{\begin{array}{c}\mathit{x}(\mathit{t})=-\frac{22+7sin(\mathit{t})}{8+4sin(\mathit{t})},\\ \mathit{y}(\mathit{t})=\frac{50+19sin(\mathit{t})}{12+6sin(\mathit{t})}\end{array}\right.$(56)

is the analytical solution (see Figure 4) to:

(57) $\left\{\begin{array}{c}tan(\mathit{t})\frac{\mathit{dx}(\mathit{t})}{\mathit{dt}}=-\frac{265}{48}-4\mathit{x}-x^2+\frac{1}{2}\mathit{y},\\ tan(\mathit{t})\frac{\mathit{dy}(\mathit{t})}{\mathit{dt}}=\frac{229}{18}-7\mathit{y}+y^2+\frac{1}{3}\mathit{x}.\end{array}\right.$(57)

Figure 4. Deformed order-1 solitary solutions $\mathit{x}(\mathit{t}),\mathit{y}(\mathit{t})$ (see eq. (56)) to (57).

7. The construction of deformed order-1 solitary solutions to (15)

7.1. The generalized operator of differentiation

The general solution to (15) is considered in the following form:

(58) $\left\{\begin{array}{c}{\hat{x}}_0:=\hat{x}(\hat{t},\mathit{c},\mathit{u},\mathit{v}),\\ {\hat{y}}_0:=\hat{y}(\hat{t},\mathit{c},\mathit{u},\mathit{v}),\end{array}\right.$(58)

where $\mathit{c}\in \mathbb{R}$ is the centre of the expansion of the solutions ${\hat{x}}_0$ and ${\hat{y}}_0$; $\mathit{u},\mathit{v}\in \mathbb{R}$ are the initial conditions $\mathit{u}={\hat{x}}_0({\hat{t}}_0,\mathit{c},\mathit{u},\mathit{v})$ and $\mathit{v}={\hat{y}}_0({\hat{t}}_0,\mathit{c},\mathit{u},\mathit{v})$.

According to (46), the simplified structure of the function $\hat{f}\left(\hat{t}\right)$ is set to:

(59) $\hat{f}\left(\hat{t}\right)=\mathit{\kappa }(\hat{t}-\mathit{\nu });$(59)

where $\mathit{\kappa }=A_1$ and $\mathit{\nu }=\frac{A_0}{A_1}$. Now, the generalized differential operator of (15) reads [25]:

(60) $D_{\mathit{cuv}}:=D_c+\frac{1}{\mathit{\kappa }(\mathit{c}-\mathit{\nu })}\left((a_0+a_1\mathit{u}+a_2{u}^2+{\epsilon }_x\mathit{v})D_u+(b_0+b_1\mathit{v}+b_2{v}^2+{\epsilon }_y\mathit{u})D_v\right)$(60)

where $D_c:=\frac{\mathrm{\partial }}{\mathrm{\partial c}}$; $D_u:=\frac{\mathrm{\partial }}{\mathrm{\partial u}}$; $D_v:=\frac{\mathrm{\partial }}{\mathrm{\partial v}}$. Then the general solution to (15) is given by:

(61) $\left\{\begin{array}{c}\widehat{x_0}(\hat{t},\mathit{c},\mathit{u},\mathit{v})={\sum }_{\mathit{j}=0}^{+\mathrm{\infty }}\frac{{(\hat{t}-\mathit{c})}^{\mathit{j}}}{\mathit{j}!}{p}_j(\mathit{c},\mathit{u},\mathit{v}),\\ \widehat{y_0}(\hat{t},\mathit{c},\mathit{u},\mathit{v})={\sum }_{\mathit{j}=0}^{+\mathrm{\infty }}\frac{{(\hat{t}-\mathit{c})}^{\mathit{j}}}{\mathit{j}!}{q}_j(\mathit{c},\mathit{u},\mathit{v}),\end{array}\right.$(61)

where

(62) $\left\{\begin{array}{c}{p}_j(\mathit{c},\mathit{u},\mathit{v}):=D_{\mathit{cuv}}^j\mathit{u},\\ q_j(\mathit{c},\mathit{u},\mathit{v}):=D_{\mathit{cuv}}^j\mathit{v}.\end{array}\right.$(62)

Solution (61) is the analytical solution to (15). However, the only solutions of interest to this study are those which can be expressed in closed form. The closed form solutions exist if the sequence of parameters $p_j(\mathit{c},\mathit{u},\mathit{v})$ and $q_j(\mathit{c},\mathit{u},\mathit{v})$ form linear recurrence sequences [8].

Lemma 7.1. The first roots of the characteristic polynomials of the linear recurrence sequences (62) are equal to 0.

Proof. The expression of the deformed order-1 solitary solutions in (27) can be rearranged as follows:

(63) $\left\{\begin{array}{c}\hat{x}\left(\hat{t}\right)={\sigma }_x\frac{\hat{t}-x_1}{\hat{t}-t_1}={\sigma }_x+\frac{{\sigma }_x(t_1-x_1)}{\hat{t}-t_1}={\sigma }_x+{\lambda }_1{\sum }_{\mathit{j}=0}^{+\mathrm{\infty }}\frac{{(\hat{t}-t_1)}^{\mathit{j}}}{\mathit{j}!}\left(\mathit{j}!{\rho }_x^j\right),\\ \hat{y}\left(\hat{t}\right)={\sigma }_y\frac{\hat{t}-y_1}{\hat{t}-t_1}={\sigma }_y+\frac{{\sigma }_y(t_1-y_1)}{\hat{t}-t_1}={\sigma }_y+{\mu }_1{\sum }_{\mathit{j}=0}^{+\mathrm{\infty }}\frac{{(\hat{t}-t_1)}^{\mathit{j}}}{\mathit{j}!}\left(\mathit{j}!{\rho }_y^j\right),\end{array}\right.$(63)

where ${\lambda }_1,{\mu }_1\in \mathbb{R}$, and the common ratios of the geometric progressions ${\rho }_x$ and ${\rho }_y$ are non-zero. Therefore, the roots of the characteristic polynomial of the linear recurrence sequence ${\sigma }_x+{\lambda }_1,{\lambda }_1{\rho }_x,{\lambda }_1{\rho }_x^2,\dots$ are $0$ and ${\rho }_x\ne 0$. Analogously, the roots of the characteristic polynomial of the linear recurrence sequence ${\sigma }_y+{\mu }_1,{\mu }_1{\rho }_y,{\mu }_1{\rho }_y^2,\dots$ are $0$ and ${\rho }_y\ne 0$.

The structure of the deformed order-1 solitary solution (27) requires that the linear recurrence sequences have an order of 2:

(64) $\left\{\begin{array}{c}{p}_j(\mathit{c},\mathit{u},\mathit{v})=\mathit{j}!\left({\lambda }_0{0}^j+{\lambda }_1{({\rho }_x(\mathit{c},\mathit{u},\mathit{v}))}^{\mathit{j}}\right),\\ q_j(\mathit{c},\mathit{u},\mathit{v})=\mathit{j}!\left({\mu }_0{0}^j+{\mu }_1{({\rho }_y(\mathit{c},\mathit{u},\mathit{v}))}^{\mathit{j}}\right),\end{array}\right.$(64)

where ${\lambda }_0,{\lambda }_1\in \mathbb{R}$; ${\rho }_x$ and ${\rho }_y$ are non-zero roots of the characteristic polynomials of the linear recurrence sequences (62).

Then, according to the basic properties of geometric progressions, the general solution to (15) reads:

(65) $\left\{\begin{array}{c}\widehat{x_0}(\hat{t},\mathit{c},\mathit{u},\mathit{v})={\lambda }_0+{\lambda }_1/(1-{\rho }_x(\mathit{c},\mathit{u},\mathit{v})(\hat{t}-\mathit{c})),\\ \widehat{y_0}(\hat{t},\mathit{c},\mathit{u},\mathit{v})={\mu }_0+{\mu }_1/(1-{\rho }_y(\mathit{c},\mathit{u},\mathit{v})(\hat{t}-\mathit{c})).\end{array}\right.$(65)

7.2. The Hankel determinants

Suppose that (64) hold true. Then, the following determinants of Hankel matrices are equal to zero at least for some values of the parameters $\mathit{c}$, $\mathit{u}$, and $\mathit{v}$:

(66) $\left\{\begin{array}{c}{H}_x(\mathit{c},\mathit{u},\mathit{v})=det\left[\begin{array}{cc}\frac{p_1}{1!}& \frac{p_2}{2!}\\ \frac{p_2}{2!}& \frac{p_3}{3!}\end{array}\right]=0;\\ H_y(\mathit{c},\mathit{u},\mathit{v})=det\left[\begin{array}{cc}\frac{q_1}{1!}& \frac{q_2}{2!}\\ \frac{q_2}{2!}& \frac{q_3}{3!}\end{array}\right]=0.\end{array}\right.$(66)

Note that higher-order Hankel determinants also vanish [25]. Furthermore, full Hankel determinants $\left[\begin{array}{ccc}\frac{p_0}{0!}& \frac{p_1}{1!}& \frac{p_2}{2!}\\ \frac{p_1}{1!}& \frac{p_2}{2!}& \frac{p_3}{3!}\\ \frac{p_2}{2!}& \frac{p_3}{3!}& \frac{p_4}{4!}\end{array}\right]$ and $\left[\begin{array}{ccc}\frac{q_0}{0!}\mathit{\hspace{0.17em}}& \frac{q_1}{1!}\mathit{\hspace{0.17em}}& \frac{q_2}{2!}\mathit{\hspace{0.17em}}\\ \frac{q_1}{1!}& \frac{q_2}{2!}& \frac{q_3}{3!}\\ \frac{q_2}{2!}& \frac{q_3}{3!}& \frac{q_4}{4!}\end{array}\right]$ are not necessary to consider since the first root ${\rho }_0$ of linear recurrence sequences (64) is zero (Lemma 7.1).

The expressions of $p_j(\mathit{c},\mathit{u},\mathit{v})$ and $q_j(\mathit{c},\mathit{u},\mathit{v})$; $\mathit{j}=1,2,3$ follow from 62. These expressions are listed in Appendix A.

Let us consider a hypothesis that the initial conditions $\mathit{u},\mathit{v}$ corresponding to order-1 solitary solutions are in the following linear relationship:

(67) $\mathit{L}:\left\{\begin{array}{c}\mathit{u}:={\alpha }_x\mathit{\omega }+{\beta }_x;\\ \mathit{v}:={\alpha }_y\mathit{\omega }+{\beta }_y,\end{array}\right.$(67)

where $\mathit{L}$ denotes a line in the phase space of initial conditions; $\mathit{\omega }$ is the parameter generating the line; ${\alpha }_x$, ${\alpha }_y$, ${\beta }_x$, ${\beta }_y\in \mathbb{R}$ are parameters defining the orientation of line $\mathit{L}$.

Note that the Hankel determinants (66) can be expressed as fifth order polynomial with respect to $\mathit{\omega }$:

(68) $\left\{\begin{array}{c}{H}_x(\mathit{c},\mathit{u},\mathit{v}{)|}_{\mathit{u},\mathit{v}\in \mathit{L}}={\tilde{H}}_{\mathit{x}}(\mathit{\omega },{\alpha }_x,{\alpha }_y,{\beta }_x,{\beta }_y)={\sum }_{\mathit{j}=0}^5{\theta }_j^{(\mathit{x})}{\omega }^j;\\ H_y(\mathit{c},\mathit{u},\mathit{v}{)|}_{\mathit{u},\mathit{v}\in \mathit{L}}={\tilde{H}}_{\mathit{y}}(\mathit{\omega },{\alpha }_x,{\alpha }_y,{\beta }_x,{\beta }_y)={\sum }_{\mathit{j}=0}^5{\theta }_j^{(\mathit{y})}{\omega }^j,\end{array}\right.$(68)

where the expressions of ${\theta }_j^{(\mathit{x})}$ and ${\theta }_j^{(\mathit{y})}$ are determined using computer algebra.

7.3. Conditions for the existence of deformed solitary solutions to 15

Coefficients ${\theta }_5^{(\mathit{x})}$ and ${\theta }_5^{(\mathit{y})}$ read:

(69) $\begin{array}{cc}{\theta }_5^{(\mathit{x})}& =4a_2{\alpha }_x^2{\alpha }_y{(a_2{\alpha }_x-b_2{\alpha }_y)}^2{\epsilon }_x;\\ {\theta }_5^{(\mathit{y})}& =4b_2{\alpha }_x{\alpha }_y^2{(a_2{\alpha }_x-b_2{\alpha }_y)}^2{\epsilon }_y.\end{array}$(69)

Equating ${\theta }_5^{(\mathit{x})}$ and ${\theta }_5^{(\mathit{y})}$ to zero yields the following relations:

(70) $\begin{array}{cc}{\alpha }_x& =\frac{b_2{\alpha }_y}{a_2};\\ {\alpha }_y& =\frac{a_2{\alpha }_x}{b_2}.\end{array}$(70)

Note that ${\theta }_5^{(\mathit{x})}$ and ${\theta }_5^{(\mathit{y})}$ may also vanish when $a_2,b_2,{\alpha }_x,{\alpha }_y,{\epsilon }_x,{\epsilon }_y$ are set to zero, however, that would result in a special case of the considered system that is either uncoupled or linear. Because of this, these cases are no longer considered.

Solution to (70) reads:

(71) ${\alpha }_x=b_2;{\alpha }_y=a_2.$(71)

Relations (71) are applied in order to simplify the coefficients ${\theta }_4^{(\mathit{x})}$ and ${\theta }_4^{(\mathit{y})}$:

(72) ${\theta }_4^{(\mathit{x})}=a_2^2{b}_2^2\left(-a_1^2{b}_2^2+4a_0{a}_2{b}_2^2-2a_2{b}_1{b}_2{\epsilon }_x+a_2^2{\epsilon }_x^2-2b_2^2{\epsilon }_x{\epsilon }_y+b_2^2{\kappa }^2\right);$(72)

(73) ${\theta }_4^{(\mathit{y})}=a_2^2{b}_2^2\left(-b_1^2{a}_2^2+4b_0{b}_2{a}_2^2-2b_2{a}_1{a}_2{\epsilon }_y+b_2^2{\epsilon }_y^2-2a_2^2{\epsilon }_y{\epsilon }_x+a_2^2{\kappa }^2\right).$(73)

Setting ${\theta }_4^{(\mathit{x})}$ to zero yields:

(74) ${\kappa }^2=a_1^2-4a_0{a}_2+\frac{{\epsilon }_x\left(2a_2{b}_1{b}_2-a_2^2{\epsilon }_x+2b_2^2{\epsilon }_y\right)}{b_2^2}.$(74)

Analogously, setting ${\theta }_4^{(\mathit{y})}$ to zero yields:

(75) ${\kappa }^2=b_1^2-4b_0{b}_2+\frac{{\epsilon }_y\left(2b_2{a}_1{a}_2-b_2^2{\epsilon }_y+2a_2^2{\epsilon }_x\right)}{a_2^2}.$(75)

Identities (74) and (75) yield the following relation:

(76) $\begin{array}{cc}\mathit{\hspace{0.17em}}& -4a_2^3{b}_2^2{a}_0+a_1^2{a}_2^2{b}_2^2-2a_1{a}_2{b}_2^3{\epsilon }_y-a_2^4{\epsilon }_x^2+2a_2^3{b}_1{b}_2{\epsilon }_x\\ \mathit{\hspace{0.17em}}& +4a_2^2{b}_0{b}_2^3-a_2^2{b}_1^2{b}_2^2+b_2^4{\epsilon }_y^2=0\end{array}$(76)

Relations (74)-(76) are applied to simplify the coefficients ${\theta }_3^{(\mathit{x})}$ and ${\theta }_3^{(\mathit{y})}$, resulting in:

(77) $\begin{array}{cc}{\theta }_3^{(\mathit{x})}& ={\epsilon }_x{\left(2a_2^2{b}_2{\beta }_x-2b_2^2{a}_2{\beta }_y+b_2{a}_2{a}_1+a_2^2{\epsilon }_x-b_2{a}_2{b}_1-{\epsilon }_y{b}_2^2\right)}^2;\\ {\theta }_3^{(\mathit{y})}& ={\epsilon }_y{\left(2a_2^2{b}_2{\beta }_x-2b_2^2{a}_2{\beta }_y+b_2{a}_2{a}_1+a_2^2{\epsilon }_x-b_2{a}_2{b}_1-{\epsilon }_y{b}_2^2\right)}^2.\end{array}$(77)

Note that due to the relation (70), coefficients ${\theta }_3^{(\mathit{x})}$ and ${\theta }_3^{(\mathit{y})}$ in (77) are proportional:

(78) ${\theta }_3^{(\mathit{y})}=\frac{{\epsilon }_y}{{\epsilon }_x}{\theta }_3^{(\mathit{x})}.$(78)

Thus, equating ${\theta }_3^{(\mathit{x})}$ and ${\theta }_3^{(\mathit{y})}$ to zero yields one equality instead of two:

(79) $2a_2{b}_2^2{\beta }_y+b_2{a}_2{b}_1+{\epsilon }_y{b}_2^2=2a_2^2{b}_2{\beta }_x+a_2{b}_2{a}_1+{\epsilon }_x{a}_2^2.$(79)

Applying all the relations mentioned above yields ${\theta }_2^{(\mathit{x})}={\theta }_2^{(\mathit{y})}={\theta }_1^{(\mathit{x})}\text{break}={\theta }_1^{(\mathit{y})}={\theta }_0^{(\mathit{x})}={\theta }_0^{(\mathit{y})}=0$ and consequently $H_x=H_y=0$.

Note that the computations presented above demonstrate that the hypothesized relation between initial conditions $\mathit{u},\mathit{v}$ (67) yields order-1 solitary solutions.

7.4. Main theorem

Derivations provided up to this point can be summarized in the following Theorem:

Theorem 7.2. The image system of Riccati equations with diffusive coupling (15) and initial conditions $\hat{x}(\mathit{c})=\mathit{u},\hat{y}(\mathit{c})=\mathit{v}$ admits the kink solitary solution (27) if:

1. Condition (76) holds true;

2. Initial conditions $\mathit{u},\mathit{v}$ lie on the following line in the phase space:

   (80) $\frac{a_2}{a_2{\beta }_x-b_2{\beta }_y}\mathit{u}+\frac{b_2}{b_2{\beta }_y-a_2{\beta }_x}\mathit{v}=1,$(80)

where ${\beta }_x,{\beta }_y$ satisfy (79).

Proof.

1. The proof is provided in Section 7.3: this condition is part of the sufficient condition for the Hankel determinants to become equal to zero.

2. Equations (67) yield the following equality:

   (81) $\frac{u-{\beta }_x}{{\alpha }_x}=\frac{v-{\beta }_y}{{\alpha }_y}.$(81)

   Basic rearrangements result in the second point of Theorem.

Let the conditions (1) and (2) of Theorem 7.2 hold true. Then, according to Lemma 7.1, the first roots of the characteristic polynomials of the linear recurrence sequences (62) are equal to 0. The remaining non-zero roots ${\rho }_x$ and ${\rho }_y$ are computed using the following characteristic equations:

(82) $\left|\begin{array}{cc}\frac{p_1}{1!}& \frac{p_2}{2!}\\ 1& {\rho }_x\end{array}\right|=0;\left|\begin{array}{cc}\frac{q_1}{1!}& \frac{q_2}{2!}\\ 1& {\rho }_y\end{array}\right|=0.$(82)

Solving (82) yields:

(83) ${\rho }_x={\rho }_y=\frac{1}{2\kappa a_2\mathit{c}}\left(a_2\left(2({\alpha }_y\mathit{\omega }+{\beta }_y)b_2-\mathit{\kappa }+b_1\right)+b_2{\epsilon }_y\right).$(83)

Parameters ${\lambda }_0$, ${\lambda }_1$, ${\mu }_0$, ${\mu }_1$ of the general solution (65) to (15) can be obtained by solving the following systems of linear equations (as seen from (64)):

(84) $\begin{array}{cc}{p}_0& =1!({\lambda }_0{0}^0+{\lambda }_1{\rho }_x^0)={\lambda }_0+{\lambda }_1;\\ p_1& =2!({\lambda }_0{0}^1+{\lambda }_1{\rho }_x^1)=2{\lambda }_1{\rho }_x;\end{array}$(84)

(85) $\begin{array}{cc}{q}_0& =1!({\mu }_0{0}^0+{\mu }_1{\rho }_y^0)={\mu }_0+{\mu }_1;\\ q_1& =2!({\mu }_0{0}^1+{\mu }_1{\rho }_y^1)=2{\mu }_1{\rho }_y;\end{array}$(85)

Note that $0^0$ is taken to be equal to $1$ [26].

Solutions to (84)–(85) are given as follows:

(86) $\begin{array}{cc}\mathit{\hspace{0.17em}}& {\lambda }_0=-\frac{a_1{b}_2+a_2{\epsilon }_x+\mathit{\kappa }{b}_2}{2a_2{b}_2};{\lambda }_1=\mathit{u}-{\lambda }_0;\\ \mathit{\hspace{0.17em}}& {\mu }_0=-\frac{b_1{a}_2+b_2{\epsilon }_y+\mathit{\kappa }{a}_2}{2b_2{a}_2};{\mu }_1=\mathit{v}-{\mu }_0.\end{array}$(86)

Using Theorem 7.2, important conclusions are made about the original system of Riccati equations with diffusive coupling (14).

Corollary 7.3. The system of Riccati equations with diffusive coupling (14) and initial conditions $\mathit{x}(t_0)=\mathit{u},\mathit{y}(t_0)=\mathit{v}$ admits the deformed kink solitary solution (11) if:

1. Condition (76) holds true;

2. Initial conditions $\mathit{u},\mathit{v}$ lie on the line 80 where ${\beta }_x,{\beta }_y$ satisfy (79);

3. The independent variable substitution $\mathit{\phi }(\mathit{t})$ reads:

   (87) $\mathit{\phi }(\mathit{t})=\mathit{\nu }+\left(\mathit{c}-\mathit{\nu }\right)exp\left(\mathit{\kappa }{\int }_{t_0}^{\mathit{t}}\frac{\mathrm{d}\tau }{\mathit{f}\left(\mathit{\tau }\right)}\right),\mathit{\phi }(t_0)=\mathit{c};$(87)

Proof. Points 1 and 2 follow directly from Theorem 7.2. To prove point 3, note that by Lemma 4.4, the independent variable substitution $\mathit{\phi }(\mathit{t})$ satisfies (16). Furthermore, $\hat{f}\left(\hat{t}\right)$ is given by (59). Combining these two equalities yields:

(88) $\frac{\mathit{d\phi }}{\mathit{dt}}=\frac{\mathit{\kappa }\left(\mathit{\phi }-\mathit{\nu }\right)}{\mathit{f}(\mathit{t})}.$(88)

Integration of the above equality with respect to $\mathit{t}$ yields (87).

The relationship between the system of Riccati equations with diffusive coupling (14), the image system (15) and their respective solutions is depicted in Figure 5.

Figure 5. Schematic summary of theorem 7.2 and corollary 7.3: the relationship between Riccati system with diffusive coupling, the image system and their kink solitary solutions.

8. The structural stability of (3)

Let us consider $\mathit{f}(\mathit{t})=1$. Then, according to Lemma 4.4 the deformation function $\mathit{\phi }(\mathit{t})$ reads:

(89) $\mathit{\phi }(\mathit{t})=exp\left(\mathit{t}\right)+\mathit{\nu }.$(89)

Inserting 89 into 11 yields:

(90) $\left\{\begin{array}{c}\mathit{x}(\mathit{t})={\sigma }_x\frac{exp\left(t\right)+\nu -x_1}{exp\left(\mathit{t}\right)+\mathit{\nu }-t_1};\\ \mathit{y}(\mathit{t})={\sigma }_y\frac{exp\left(t\right)+\nu -y_1}{exp\left(\mathit{t}\right)+\mathit{\nu }-t_1},\end{array}\right.$(90)

which are the non-deformed order–1 solitary solutions 9 to 8 (at $\mathit{\eta }=1$ and $t_0=0$). In other words, Theorem 7.2 implies that Hypothesis 3.1 does not hold true.

This result is completely counter-intuitive. The transition from the hepatitis C virus infection model (3) to the model of two Riccati equations coupled with diffusive terms (8) can be executed only by erasing the multiplicative terms in (3). However, the necessary conditions for the existence of non-deformed order–1 solitary solutions to (3) require the identities (10) to hold true. But identities (10) eliminate the quadratic terms from (8). In other words, Riccati equations in (8) are reduced to linear differential equations. Therefore, the limiting transition in the hepatitis C virus infection model does not allow the existence of non-deformed order–1 solitary solutions to (8). However, Theorem 7.2 yields an opposite result. It is clear that non-deformed order-1 solitary solutions to (8) indeed exist.

Structural stability is a fundamental property of a dynamical system which means that the qualitative behaviour of the trajectories is unaffected by small perturbations [27]. The limiting transition in (3) accompanied by the necessary and sufficient conditions for the existence of non-deformed order-1 solitary solutions does not produce correct results. In other words, the perturbation from (8) to (3) changes the existence of solitary solutions if executed by adding the multiplicative coupling terms (even with infinitesimal, but non-zero weight coefficients).

Corollary 8.1. The hepatitis C virus infection model 1 is structurally unstable with respect to the multiplicative coupling terms.

9. Concluding remarks

It is a common feature of nonlinear dynamical systems that a solitary solution form a separatrix in the space of system parameters and initial conditions [6]. Therefore, the existence of solitary solutions may help to understand the global dynamics of the system [28]. In particular, kink solitary solutions to the hepatitis C virus infection model can be either in a linear or in a hyperbolic relationship ([7]). Thus, a large perturbation in the population of hepatitis infected cells does not necessarily lead to a large change in uninfected cells (if the relationship is hyperbolic) ([7]).

However, the results of this paper imply that the situation with the order-1 solitary solutions to the hepatitis C virus infection model is even more complex. The limit transition from the hepatitis C virus infection model to the model of Riccati equations coupled with the diffusive terms does not destroy the existence of the order-1 solitary solutions. Furthermore, this controversy does not occur in the other limit transition from the hepatitis C virus infection model to the model of Riccati equations coupled with the multiplicative terms (setting $a_4$ and $b_4$ to zero in (3) does not destroy the structure of Riccati equations).

Apart from the biological relevance, this paper provides a significant contribution to the mathematical theory of solitary solutions to nonlinear differential equations. It appears that the role of diffusive and multiplicative coupling terms in the dynamics of nonlinear systems is very different not only from the point of view for the existence and the stability of equilibria, but also with respect to the existence of solitary solutions to the coupled equations. Further extension of the complexity of the coupled model to additional spatial dimensions (for example, coupled systems on heterogeneous graphs) and towards more complex nodal nonlinearity (beyond Riccati nomenclature) remains a definite objective of future research.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Appendix A.

Coefficients ${\mathit{p}}_1,{\mathit{p}}_2,{\mathit{p}}_3$ and ${\mathit{q}}_1,{\mathit{q}}_2,{\mathit{q}}_3$

(A1) $p_1=\frac{1}{\mathit{\kappa }\left(\mathit{c}-\mathit{\nu }\right)}\left(a_0+a_1\mathit{u}+a_2{u}^2+{\epsilon }_x\mathit{v}\right);$(A1)

(A2) $q_1=\frac{1}{\mathit{\kappa }\left(\mathit{c}-\mathit{\nu }\right)}\left(b_0+b_1\mathit{v}+b_2{v}^2+{\epsilon }_y\mathit{u}\right);$(A2)

(A3) $\begin{array}{cc}{p}_2& =\frac{1}{{\kappa }^2{\left(\mathit{c}-\mathit{\nu }\right)}^2}(2u^3{a_2}^2-a_2\left(\mathit{\kappa }-3a_1\right)u^2\\ \mathit{\hspace{0.17em}}& +\left(\left(2\mathit{v}{a}_2+{\epsilon }_y\right){\epsilon }_x-\mathit{\kappa }{a}_1+2a_0{a}_2+{a_1}^2\right)\mathit{u}\\ \mathit{\hspace{0.17em}}& +\left(\mathit{v}{a}_1+v^2{b}_2+\left(-\mathit{\kappa }+b_1\right)\mathit{v}+b_0\right){\epsilon }_x-a_0\left(\mathit{\kappa }-a_1\right));\end{array}$(A3)

(A4) $\begin{array}{cc}{q}_2& =\frac{1}{{\kappa }^2{\left(\mathit{c}-\mathit{\nu }\right)}^2}(2v^3{b_2}^2-b_2\left(\mathit{\kappa }-3b_1\right)v^2\\ \mathit{\hspace{0.17em}}& +\left(\left(2\mathit{v}{b}_2+{\epsilon }_x\right){\epsilon }_y-\mathit{\kappa }{b}_1+2b_0{b}_2+{b_1}^2\right)\mathit{v}\\ \mathit{\hspace{0.17em}}& +\left(\mathit{v}{b}_1+u^2{a}_2+\left(-\mathit{\kappa }+a_1\right)\mathit{u}+a_0\right){\epsilon }_y-b_0\left(\mathit{\kappa }-b_1\right));\end{array}$(A4)

(A5) $\begin{array}{cc}{p}_3& =\frac{1}{{\kappa }^3{\left(\mathit{c}-\mathit{\nu }\right)}^3}(\mathit{v}\left(2\mathit{v}{a}_2+{\epsilon }_y\right){{\epsilon }_x}^2+\left(\left(8\mathit{v}{a_2}^2+3a_2{\epsilon }_y\right)u^2\right.\\ \mathit{\hspace{0.17em}}& +\left(2b_2{v}^2{a}_2+\left(\left(2b_1-6\mathit{\kappa }+8a_1\right)a_2+2b_2{\epsilon }_y\right)\mathit{v}+2b_0{a}_2\right.\\ \mathit{\hspace{0.17em}}& \left.-3{\epsilon }_y\left(\mathit{\kappa }-\frac{2}{3}{a}_1-\frac{1}{3}{b}_1\right)\right)\mathit{u}+2{b_2}^2{v}^3-3\left(\mathit{\kappa }-\frac{1}{3}{a}_1-b_1\right)b_2{v}^2\\ \mathit{\hspace{0.17em}}& +\left(4a_0{a}_2+{a_1}^2+\left(b_1-3\mathit{\kappa }\right)a_1+{b_1}^2-3b_1\mathit{\kappa }+2b_2{b}_0+2{\kappa }^2\right)\mathit{v}\\ \mathit{\hspace{0.17em}}& \left.+{\epsilon }_y{a}_0-3b_0\mathit{\kappa }+b_0{a}_1+b_1{b}_0\right){\epsilon }_x+6\left(a_2{u}^2+a_1\mathit{u}+a_0\right)\\ \mathit{\hspace{0.17em}}& \times \left(u^2{a_2}^2-a_2\left(\mathit{\kappa }-a_1\right)\mathit{u}+\frac{1}{3}{a}_0{a}_2+\frac{1}{3}\left(\mathit{\kappa }-a_1\right)\left(\mathit{\kappa }-\frac{1}{2}{a}_1\right)\right));\end{array}$(A5)

(A6) $\begin{array}{cc}{q}_3& =\frac{1}{{\kappa }^3{\left(\mathit{c}-\mathit{\nu }\right)}^3}(\mathit{u}\left(2\mathit{v}{b}_2+{\epsilon }_x\right){{\epsilon }_y}^2+\left(\left(8\mathit{v}{b_2}^2+3b_2{\epsilon }_x\right)v^2\right.\\ \mathit{\hspace{0.17em}}& +\left(2a_2{u}^2{b}_2+\left(\left(2a_1-6\mathit{\kappa }+8b_1\right)b_2+2a_2{\epsilon }_x\right)\mathit{u}+2a_0{b}_2\right.\\ \mathit{\hspace{0.17em}}& \left.-3{\epsilon }_x\left(\mathit{\kappa }-\frac{2}{3}{b}_1-\frac{1}{3}{a}_1\right)\right)\mathit{v}+2{a_2}^2{u}^3-3\left(\mathit{\kappa }-\frac{1}{3}{b}_1-a_1\right)a_2{u}^2\\ \mathit{\hspace{0.17em}}& +\left(4b_0{b}_2+{b_1}^2+\left(a_1-3\mathit{\kappa }\right)b_1+{a_1}^2-3a_1\mathit{\kappa }+2a_2{a}_0+2{\kappa }^2\right)\mathit{u}\\ \mathit{\hspace{0.17em}}& \left.+{\epsilon }_x{b}_0-3a_0\mathit{\kappa }+a_0{b}_1+a_1{a}_0\right){\epsilon }_y+6\left(b_2{v}^2+b_1\mathit{v}+b_0\right)\\ \mathit{\hspace{0.17em}}& \times \left(v^2{b_2}^2-b_2\left(\mathit{\kappa }-b_1\right)\mathit{v}+\frac{1}{3}{b}_0{b}_2+\frac{1}{3}\left(\mathit{\kappa }-b_1\right)\left(\mathit{\kappa }-\frac{1}{2}{b}_1\right)\right)).\end{array}$(A6)