7.2
CiteScore
3.7
Impact Factor
Generic selectors
Exact matches only
Search in title
Search in content
Post Type Selectors
Search in posts
Search in pages
Filter by Categories
ABUNDANCE ESTIMATION IN AN ARID ENVIRONMENT
Case Study
Correspondence
Corrigendum
Editorial
Full Length Article
Invited review
Letter to the Editor
Original Article
Retraction notice
REVIEW
Review Article
SHORT COMMUNICATION
Short review
7.2
CiteScore
3.7
Impact Factor
Generic selectors
Exact matches only
Search in title
Search in content
Post Type Selectors
Search in posts
Search in pages
Filter by Categories
ABUNDANCE ESTIMATION IN AN ARID ENVIRONMENT
Case Study
Correspondence
Corrigendum
Editorial
Full Length Article
Invited review
Letter to the Editor
Original Article
Retraction notice
REVIEW
Review Article
SHORT COMMUNICATION
Short review
View/Download PDF

Translate this page into:

Original article
07 2023
:35;
102682
doi:
10.1016/j.jksus.2023.102682

Solitary solutions to a metastasis model represented by two systems of coupled Riccati equations

Center for Nonlinear Systems, Kaunas University of Technology, Studentu 50-147, Kaunas LT-51368, Lithuania
Department of Software Engineering, Kaunas University of Technology, Studentu 50-415, Kaunas LT-51368, Lithuania
Urology Clinic, Lithuanian University of Health Sciences, Eiveniu 2, Kaunas LT-50009, Lithuania

⁎Corresponding author. inga.timofejeva@ktu.edu (I. Timofejeva)

Disclaimer:
This article was originally published by Elsevier and was migrated to Scientific Scholar after the change of Publisher.

Peer review under responsibility of King Saud University.

Abstract

Solitary solutions to a two-tumor metastasis model represented by a system of multiplicatively coupled Riccati equations are considered in this paper. The interaction between tumors is modelled by a one-sided diffusive coupling when one coupled Riccati system influences the other, but not the opposite. Necessary and sufficient conditions for the existence of solitary solutions to the composite system of Riccati equations are derived in the explicit form. Computational experiments are used to demonstrate the transitions from one steady-state to another via non-monotonous trajectories.

Keywords

Nonlinear differential equation
Analytical solution
Tumor model
35C08
34A25
92C50
PubMed
1

1 Introduction

Mathematical modeling of biomedical processes, and cancer, in particular, has come to the forefront of research in recent decades (Altrock et al., 2015). A conceptual framework for cancer biology aims to identify the structure, the dynamics, and the rules governing the regulation of the behavior of the system (Medina, 2018).

A vast plethora of cancer models have been introduced into the scientific literature during the last years (Magi et al., 2017). In general, these models can be classified according to the mathematical modeling approach: the top-down approach and the bottom-up approach (Shahzad and Loor, 2012). Also, mathematical models of cancer are built at different scales with different degrees of detail (Medina, 2018): very detailed models of metabolic pathways, middle scale models based on fluxomics, and global phenomenological scale cancer models based on the interaction of populations of cells.

Global phenomenological scale cancer modeling approaches can be characterized either as continuous models or discrete models (with various subclasses under each type) (Simmons et al., 2017). Continuum models generate solutions as concentrations of different types of cells. Usually, they comprise a system of differential equations (though agent-based, multi-scale, lattice-based, image-based models are widely used also) (Simmons et al., 2017). Ordinary differential equations are employed if the variation of cell populations is investigated in time. Partial differential equations are employed to describe the variation both in time and in space.

Mathematical models of a single cancer tumor comprising a system of two ordinary differential equations coupled by multiplicative terms are widely used for the phenomenological description of the interaction between healthy and cancer cells (Liu and Yang, 2016). Such types of models comprising two Riccati type nonlinear coupled differential equations are used for the description of prostate cancer treatment with androgen deprivation therapy (Zazoua and Wang, 2019), cancer stem cell-targeted immunotherapy (Sigal et al., 2019; Alqudah, 2020), the maximization of viability time in general cancer therapy (Bratus et al., 2017).

The phenomenological model of such systems in mathematical terms reads:

(1)
x = a 0 + a 1 x + a 2 x 2 + λ x xy ; y = b 0 + b 1 y + b 2 y 2 + λ y xy , where a 0 , a 1 , a 2 0 , b 0 , b 1 , b 2 0 R and λ x , λ y - coefficients of multiplicative terms.

A formal classification of multiple cancers in humans is first done in Moertel (1977). According to this classification, multiple primary cancers of the first type occur in organs with the same histology. The second type of multiple primary cancers originate from different tissues or organs. The third type of multiple cancers comprise three or more cancers at different tissues and organs that concurrently coexist with the first type of cancers (Moertel, 1977). Cancers of the first type are further subdivided into three groups: 1A – two or more cancers that occur in the same tissue or organ; 1B – cancers in a joint tissue shared by different organs; 1C – cancers occurring in bilaterally paired organs (Moertel, 1977).

Multiple primary cancers (MPC) are classified as synchronous and metachronous cancers (Jena et al., 2016). Among patients with MPC, triple cancers occur in 0.5 %, and quadruple or quintuple cancers in less than 0.1 % of the population (Kim et al., 2010; Kim et al., 2013). For metachronous gastric cancer patients, the most common cancer sites are: colorectum, thyroid, lung, kidney, and breast (in descending order) (Kim et al., 2013). For synchronous gastric cancer, a different order is observed: head and neck, esophagus, lung, and kidney (Kim et al., 2013). Females are at a higher risk of being diagnosed with a first primary cancer of breast or thyroid cancer at a younger age – but those who survived their first primary cancer are at a higher risk of developing a metachronous second primary cancer after that (Jena et al., 2016).

The life expectancy for a person diagnosed with multiple cancer is significantly shorter compared to a person only with a single isolated cancer (Jena et al., 2016). For instance, MPC developed in 10.4% of patients during five years after stomach cancer treatment (Jena et al., 2016). The probability of survival of such patients with MPC was considerably less than of patients with only a single stomach cancer (Jena et al., 2016). The majority of patients with synchronous MPC died during two years after the removal of the stomach. The survival rate of patients with metachronous MPC was also low – but many of those patients did survive for five years. The percentages of patients who did survive for five years after the removal of the stomach where: 76.5% for patients without MPC; 67.5% for patients with metachronous MPC; 34.1% for patients with synchronous MPC (Jena et al., 2016).

A mathematical model for the immune-mediated theory of metastasis is recently presented in Rhodes and Hillen (2019). The interaction between different tumors is modelled by using diffusive terms representing the migration of cancer cells between primary and secondary tumors in Rhodes and Hillen (2019). This corresponds well to the mathematical model of migration in diffusively coupled prey-predator system on heterogenous oriented graphs (Nagatani, 2020).

The main objective of this paper is to explore the existence of solitary solutions in a paradigmatic cancer metastasis model where the model of a single tumor is represented by a system of coupled ordinary differential equations with multiplicative terms, and the interaction between primary and secondary tumors is represented by a diffusive term.

In mathematical terms, this model can be represented by the following system of nonlinear differential equations:

(2)
x = a 0 + a 1 x + a 2 x 2 + λ x xy ; y = b 0 + b 1 y + b 2 y 2 + λ y xy + γ w - y ; z = c 0 + c 1 z + c 2 z 2 + λ z zw ; w = d 0 + d 1 w + d 2 w 2 + λ w zw .

This model uses a system of coupled Riccati type differential equations to characterize the dynamics of a single tumor (Liu and Yang, 2016; Zazoua and Wang, 2019; Sigal et al., 2019; Alqudah, 2020; Bratus et al., 2017). On the other hand, the diffusive terms coupling different tumors fall in line with the model of metastasis presented in Rhodes and Hillen (2019). Moreover, Eq. (2) does generalize the model of diffusively coupled prey-predator system on heterogenous graphs (Nagatani, 2020) because the governing differential equations representing the dynamics of healthy and cancer cells in each tumor are nonlinear Riccati differential equations (linear differential equations are used in Nagatani (2020)). The model presented by Eq. (2) also can be interpreted as the phenomenological generalization of the Master–Slave model comprising two discrete logistic maps (Ramirez et al., 2018). The dissipative coupling between different tumors is a one-way coupling. In other words, the primary (the “Master”) tumor does have a direct impact on the growth of the secondary (the “Slave”) tumor, but not vice versa. As mentioned previously, such dissipative connection between tumors does represent the model of the metastasis discussed in Rhodes and Hillen (2019).

Such kind of generalization of the metastasis model on heterogenous graph poses an important theoretical challenge from the mathematical point of view. The existence of solitary solutions in an androgen-deprivation prostate cancer treatment model (Telksnys et al., 2020), a hepatitis C evolution model (Telksnys et al., 2018), or a model of a biological neuron (Telksnys et al., 2019) always allows to construct the global view of the system dynamics in the phase space of system parameters and initial conditions. The matter of fact is, that solitary solutions to these models do usually represent separatrix between different basin boundaries (Telksnys et al., 2019; Telksnys et al., 2018; Telksnys et al., 2020). Therefore, the ability to identify solitary solutions (and the corresponding basin boundaries) would provide an important insight into the dynamics of the paradigmatic cancer metastasis model.

2

2 Preliminaries

2.1

2.1 The concept of solitary solutions

Analytical functions of the form:

(3)
x ( t ) = σ k = 1 n ( exp ( η ( t - t 0 ) ) - x k ) k = 1 n ( exp ( η ( t - t 0 ) ) - t k ) , are considered solitary functions in this paper; here n N represents the order of the solitary function; t 0 , σ , η R are real constants; x k , t k C parameters that depend only on the initial conditions.

The substitution:

(4)
t ^ = exp ( η ( t - t 0 ) ) , is used to rewrite (3) in a more concise form:
(5)
x ( t ) = x ln t ^ η + t 0 = x ^ t ^ = σ X t ^ T t ^ ,
where functions X , T read:
(6)
X ( θ ) = k = 1 n θ - x k ; T ( θ ) = k = 1 n θ - t k .

2.2

2.2 The inverse balancing technique

The primary question that has to be answered when considering solitary solutions to a particular system is whether they exist as solutions, and if so, what conditions must the system parameters and solution parameters satisfy. The inverse balancing technique, presented in detail in Navickas et al. (2016), given a system or differential equations, both determines the maximum solitary solution order and derives the necessary and sufficient conditions for their existence.

The fundamental idea of the technique is to insert the expression (5) as the ansatz into the system of differential equations. After simplifying, a linear system of equations with respect to the equation parameters is obtained. The conditions under which this system is non-singular coincide with the necessary existence solutions for the solitary solutions.

This approach differs significantly from the direct construction: in the latter, after inserting the ansatz into the system of differential equations, the resulting nonlinear equations are solved for the solution, rather than the system parameters. This yields a high-order system of equations that cannot always be solved, especially if all solutions are to be obtained. Due to these complications, we forgo the direct approach and use the inverse balancing technique to analyze (2).

3

3 Existence of solitary solutions in the Riccati system

The substitution (4) can be used on (2), resulting in:

(7)
η t ^ x ^ t ^ = a 0 + a 1 x ^ + a 2 x ^ 2 + λ x x ^ y ^ ;
(8)
η t ^ y ^ t ^ = b 0 + b 1 y ^ + b 2 y ^ 2 + λ y x ^ y ^ + γ w ^ - y ^ ;
(9)
η t ^ z ^ t ^ = c 0 + c 1 z ^ + c 2 z ^ 2 + λ z z ^ w ^ ;
(10)
η t ^ w ^ t ^ = d 0 + d 1 w ^ + d 2 w ^ 2 + λ w z ^ w ^ .

For the derivations presented, the following functions have to be defined:

(11)
x t = x ^ t ^ = σ 1 X T ; y t = y ^ t ^ = σ 2 Y T ;
(12)
z t = z ^ t ^ = σ 3 Z T ; w t = w ^ t ^ = σ 4 W T ;
(13)
X = X t ^ = t ^ - x 1 t ^ - x 2 t ^ - x n = t ^ n + χ n - 1 t ^ n - 1 + + χ 0 ;
(14)
Y = Y t ^ = t ^ - y 1 t ^ - y 2 t ^ - y n = t ^ n + θ n - 1 t ^ n - 1 + + θ 0 ;
(15)
Z = Z t ^ = t ^ - z 1 t ^ - z 2 t ^ - z n = t ^ n + κ n - 1 t ^ n - 1 + + κ 0 ;
(16)
W = W t ^ = t ^ - w 1 t ^ - w 2 t ^ - w n = t ^ n + ρ n - 1 t ^ n - 1 + + ρ 0 ;
(17)
T = T t ^ = t ^ - t 1 t ^ - t 2 t ^ - t n = t ^ n + æ n - 1 t ^ n - 1 + + æ 0 ;
(18)
X = X t ^ t ^ ; Y = Y t ^ t ^ ;
(19)
Z = Z t ^ t ^ ; W = W t ^ t ^ ; T = T t ^ t ^ .
Coefficients æ k , χ k , θ k , κ k , ρ k , k = 1 , , n - 1 are, in general, complex constants. The order n of polynomials X , Y , Z , W , T coincides with the order of the solutions x , y , z , w .

Inserting the ansatz (11)–(12) into (7)–(10) yields:

(20)
η t ^ σ 1 X T - XT T 2 = a 0 + a 1 σ 1 X T + a 2 σ 1 2 X 2 T 2 + λ x σ 1 σ 2 XY T 2 ;
(21)
η t ^ σ 2 Y T - YT T 2 = b 0 + b 1 - γ σ 2 Y T + b 2 σ 2 2 Y 2 T 2 + λ y σ 1 σ 2 XY T 2 + γ σ 4 W T ;
(22)
η t ^ σ 3 Z T - ZT T 2 = c 0 + c 1 σ 3 Z T + c 2 σ 3 2 Z 2 T 2 + λ z σ 3 σ 4 ZW T 2 ;
(23)
η t ^ σ 4 W T - WT T 2 = d 0 + d 1 σ 4 W T + d 2 σ 4 2 W 2 T 2 + λ w σ 3 σ 4 ZW T 2 .

The first two equations can be multiplied by T 2 X and T 2 Y , respectively; the last two by T, resulting in a simplified version of (20)–(23):

(24)
η t ^ σ 1 X T X - a 0 T 2 X = η t ^ σ 1 T + a 1 σ 1 T + a 2 σ 1 2 X + λ x σ 1 σ 2 Y ;
(25)
- η t ^ σ 2 YT T - b 2 σ 2 2 Y 2 T - λ y σ 1 σ 2 XY T = - η t ^ σ 2 Y + b 0 T + b 1 - γ σ 2 Y + γ σ 4 W ;
(26)
η t ^ σ 3 Z T Z - c 0 T 2 Z = η t ^ σ 3 T + c 1 σ 3 T + c 2 σ 3 2 Z + λ z σ 3 σ 4 W ;
(27)
- η t ^ σ 4 WT T - d 2 σ 4 2 W 2 T - λ w σ 3 σ 4 ZW T = - η t ^ σ 4 W + d 0 T + d 1 σ 4 W .

Each of the above equations can only hold if both right and left terms are polynomials of order n, which means the left side denominators must cancel out in each case:

(28)
η t ^ σ 1 X - a 0 T = σ 1 α 1 X ;
(29)
- η t ^ σ 2 T - b 2 σ 2 2 Y - λ y σ 1 σ 2 X = σ 2 β 1 T ,
(30)
η t ^ σ 3 Z - c 0 T = σ 3 α 2 Z ;
(31)
- η t ^ σ 4 T - d 2 σ 4 2 W - λ w σ 3 σ 4 Z = σ 4 β 2 T ,
where α 1 , α 2 , β 1 , β 2 R , α 1 , α 2 , β 1 , β 2 0 are constants.

The above relations coincide with the necessary existence conditions for solitary solutions to (2). They can be further simplified as:

(32)
a 0 T + α 1 σ 1 X - η t ^ σ 1 X = 0 ;
(33)
η t ^ T + β 1 T + λ y σ 1 X + b 2 σ 2 Y = 0 ,
(34)
c 0 T + α 2 σ 3 Z - η t ^ σ 3 Z = 0 ;
(35)
η t ^ T + β 2 T + λ w σ 3 Z + d 2 σ 4 W = 0 .

Combining (28)–(31) with (24)–(27) results in:

(36)
α 1 T = η t ^ T + a 1 T + a 2 σ 1 X + λ x σ 2 Y ;
(37)
σ 2 β 1 Y = - η t ^ σ 2 Y + b 0 T + b 1 - γ σ 2 Y + γ σ 4 W ;
(38)
α 2 T = η t ^ T + c 1 T + c 2 σ 3 Z + λ z σ 4 W ;
(39)
σ 4 β 2 W = - η t ^ σ 4 W + d 0 T + d 1 σ 4 W .

Obtained Eqs. (36)–(39) are sufficient conditions for solutions (11) to exist in the system (2). These conditions can be further rearranged:

(40)
η t ^ T + ( a 1 - α 1 ) T + a 2 σ 1 X + λ x σ 2 Y = 0 ;
(41)
γ σ 4 W + ( b 1 - γ - β 1 ) σ 2 Y + b 0 T - η σ 2 t ^ Y = 0 ;
(42)
η t ^ T + ( c 1 - α 2 ) T + c 2 σ 3 Z + λ z σ 4 W = 0 ;
(43)
( d 1 - β 2 ) σ 4 W + d 0 T - η σ 4 t ^ W = 0 .
Lemma 3.1

Solutions of the form (11)–(12) exist in the system (2) iff (32)–(35) and (40)–(43) hold true.

4

4 Maximum possible order of solitary solutions to the Riccati system

The inverse balancing technique, discussed in detail in Section 2.2, is used to determine the maximum possible order solitary solutions to the considered Riccati system (2).

The following isomorphisms are introduced:

(44)
T T = 1 , æ n - 1 , , æ 0 ;
(45)
X X l = 1 , χ n - 1 , , χ 0 ;
(46)
Y Y l = 1 , θ n - 1 , , θ 0 ;
(47)
Z Y l = 1 , κ n - 1 , , κ 0 ;
(48)
W Y l = 1 , ρ n - 1 , , ρ 0 ;
(49)
t ^ T t ^ T = n , ( n - 1 ) æ n - 1 , , æ 1 , 0 ;
(50)
t ^ X t ^ X = n , ( n - 1 ) χ n - 1 , , χ 1 , 0 ;
(51)
t ^ Y t ^ Y = n , ( n - 1 ) θ n - 1 , , θ 1 , 0 ,
(52)
t ^ Z t ^ Z = n , ( n - 1 ) κ n - 1 , , κ 1 , 0 ,
(53)
t ^ W t ^ W = n , ( n - 1 ) ρ n - 1 , , ρ 1 , 0 ,
the functions X , Y , Z , W , T have been noted in (13)–(17). Next, the following auxiliary parameters are introduced in order to simplify the further derivations: A 1 l , A 2 l , L 1 l , L 2 l , L 3 l , N 1 l , N 2 l , N 3 l , K 1 l , K 2 l , K 3 l R ; l = 1 , 2 .

Applying the above identities, previously derived conditions (32)–(35) and (40)–(43) can be restated via the defined vectors:

(54)
T + A 11 X + A 21 t ^ X = 0 ; T + A 12 Z + A 22 t ^ Z = 0 ;
(55)
t ^ T + L 11 T + L 21 X + L 31 Y = 0 ; t ^ T + L 12 T + L 22 Z + L 32 W = 0 ;
(56)
t ^ T + N 11 T + N 21 X + N 31 Y = 0 ; t ^ T + N 12 T + N 22 Z + N 32 W = 0 ;
(57)
W + K 11 Y + K 21 T + K 31 t ^ Y = 0 ; K 12 W + K 22 T + K 32 t ^ W = 0 .
Parameters used in the above equations read:
(58)
( i ) A kl : A 11 = α 1 σ 1 a 0 , A 12 = α 2 σ 3 c 0 ; A 21 = - η σ 1 a 0 , A 22 = - η σ 3 c 0 ; ( ii ) L ml : L 11 = β 1 η , L 12 = β 2 η , L 21 = λ y σ 1 η ; L 22 = λ w σ 3 η , L 31 = b 2 σ 2 η , L 32 = d 2 σ 4 η ; ( iii ) N ml : N 11 = a 1 - α 1 η , N 12 = c 1 - α 2 η , N 21 = a 2 σ 1 η ; N 22 = c 2 σ 3 η , N 31 = λ x σ 2 η , N 32 = λ z σ 4 η ; ( iv ) K ml : K 11 = ( b 1 - γ - β 1 ) σ 2 γ σ 4 , K 12 = ( d 1 - β 2 ) σ 4 , K 21 = b 0 γ σ 4 ; K 22 = d 0 , K 31 = - η σ 2 γ σ 4 , K 32 = - η σ 4 ,
where m = 1 , 2 , 3 ; l , k = 1 , 2 .

All of the parameters defined above can be placed in one of three exclusive groups which can be represented by a schematic diagram in Fig. 1:

  • a k , b k , c k , d k λ x , λ y , λ z , λ w and γ , k = 0 , 1 , 2 are parameters of the Riccati system (2).

  • η , t 0 , σ k , æ m , χ m , θ m , κ m and ρ m where m = 1 , , n - 1 ; k = 1 , 2 , 3 , 4 are parameters of the solitary solutions.

  • α l , β l , A kl , L ml , N ml and K ml , where m = 1 , 2 , 3 ; l , k = 1 , 2 are auxiliary parameters that neither belong to the solution nor the system of differential equations. The first two, α l and β l , are applied in the derivation of conditions (54)–(57) (specifically (28)–(31)). The remaining parameters A kl , L ml , N ml , K ml relate the Riccati system and its solitary solution parameters via necessary and sufficient conditions (54)–(57).

A schematic diagram portraying the relationship between different groups of parameters.
Fig. 1
A schematic diagram portraying the relationship between different groups of parameters.

The solitary solution and auxiliary parameters can be derived as a solution to the algebraic Eqs. (54)–(57). The Riccati system parameters are then obtained as:

(59)
a 0 = - η σ 1 A 21 , c 0 = - η σ 3 A 22 , a 1 = η N 11 - A 11 A 21 , c 1 = η N 12 - A 12 A 22 , a 2 = N 21 η σ 1 , c 2 = N 22 η σ 3 , λ x = N 31 η σ 2 , λ z = N 32 η σ 4 , b 0 = - η σ 2 K 21 K 31 , d 0 = K 22 , b 1 = - η K 11 K 31 + σ 2 K 31 σ 4 + L 11 , d 1 = K 12 σ 4 + η L 12 , b 2 = L 31 η σ 2 , d 2 = L 32 η σ 4 , λ y = L 21 η σ 1 , λ w = L 22 η σ 3 , γ = - η σ 2 K 31 σ 4 . Using the vector definitions (44), (45), (47), (50) and (52) and taking each element of (54) yields a set of algebraic equations with unknowns A 1 l , A 2 l , æ k , χ k and κ k , k = 0 , , n - 1 , l = 1 , 2 :
(60)
1 + A 11 + nA 21 = 0 ; 1 + A 12 + nA 22 = 0 ; æ k + A 11 χ k + A 21 k χ k = 0 , k = n - 1 , , 0 , æ k + A 12 κ k + A 22 k κ k = 0 , k = n - 1 , , 0 .
The solution of the above system reads:
(61)
A 2 l = - 1 + A 1 l n , χ k = æ k 1 + ( n - k ) A 21 = p k 1 æ k ,
(62)
κ k = æ k 1 + ( n - k ) A 22 = p k 2 æ k , A 1 l , æ k R .
Note that p kl is defined as p kl = 1 1 + ( n - k ) A 2 l , k = 0 , , n - 1 , l = 1 , 2 in the above.

Similarly to (60), using (44)–(49) in conjunction with (55) yields the following set of linear equations with L 1 l , L 2 l , L 3 l as unknowns:

(63)
n + L 11 + L 21 + L 31 = 0 ; n + L 12 + L 22 + L 32 = 0 ; k æ k + L 11 æ k + L 21 χ k + L 31 θ k = 0 , k = n - 1 , , 0 , k æ k + L 12 æ k + L 22 κ k + L 32 ρ k = 0 , k = n - 1 , , 0 .

Let:

(64)
θ k = h k 1 æ k , ρ k = h k 2 æ k , h kl C , k = 0 , , n - 1 , l = 1 , 2 .

Note that if æ k = 0 , the last two equations of (63) become identities. Otherwise, cancelling æ k 0 from the last two equations of (63) results in:

(65)
n + L 11 + L 21 + L 31 = 0 ; n + L 12 + L 22 + L 32 = 0 ; k + L 11 + L 21 p k 1 + L 31 h k 1 = 0 , k = n - 1 , , 0 , k + L 12 + L 22 p k 2 + L 32 h k 2 = 0 , k = n - 1 , , 0 .

Applying (44)–(49)(56) yields another set of linear equations with unknowns N 1 l , N 2 l , N 3 l :

(66)
n + N 11 + N 21 + N 31 = 0 ; n + N 12 + N 22 + N 32 = 0 ; k + N 11 + N 21 p k 1 + N 31 h k 1 = 0 , k = n - 1 , , 0 , k + N 12 + N 22 p k 2 + N 32 h k 2 = 0 , k = n - 1 , , 0 . Note that the two above systems (65) and (66) coincide fully, thus their solutions will also be derived in an analogous manner.

Finally, applying (44), (46), (48), (51) and (53)–(57) results in a set of linear equations with respect to K 1 l , K 2 l and K 3 l , l = 1 , 2 :

(67)
1 + K 11 + K 21 + nK 31 = 0 ; K 12 + K 22 + nK 32 = 0 ; h k 2 + K 11 h k 1 + K 21 + K 31 kh k 1 = 0 , k = n - 1 , , 0 , K 1 l h k 2 + K 22 + K 32 kh k 2 = 0 , k = n - 1 , , 0 .

4.1

4.1 First order solitary solutions to (2)

First order solitary solutions (kink solutions) will be considered in this subsection. The selection n = 1 transforms systems (65), (66) into a singular form; the solutions in terms of L 3 l , N 3 l read:

(68)
L 1 l = ( 1 + L 3 l ) p 0 l - L 3 l h 0 l 1 - p 0 l , L 2 l = - ( 1 + L 3 l ) + L 3 l h 0 l 1 - p 0 l , L 3 l R ;
(69)
N 1 l = ( 1 + N 3 l ) p 0 l - N 3 l h 0 l 1 - p 0 l , N 2 l = - ( 1 + N 3 l ) + N 3 l h 0 l 1 - p 0 l , N 3 l R ,
with l = 1 , 2 . While the expressions of L kl , N kl coincide, in general L kl N kl , since the values L 3 l and N 3 l can be selected arbitrarily and do not depend on each other.

The system (67) yields the expressions for K 11 , K 12 , K 21 , K 12 :

(70)
K 11 = h 02 - nK 31 - 1 1 - h 01 , K 12 = - nK 32 1 - h 02 , K 21 = ( nK 31 + 1 ) h 01 - h 02 1 - h 01 , K 22 = nK 32 h 02 1 - h 02 , where K 31 , K 32 R .
Corollary 4.1

The Riccati system (2) admits the first order solitary solutions (kink solutions, n = 1 ) without any constraints on Riccati system parameters or solitary solution parameters.

4.2

4.2 Higher order solitary solutions to (2)

Higher order solitary solutions are obtained by selecting n 2 . In this case (65), (66) are no longer singular and can only admit one solution:

(71)
L 1 l = N 1 l = ( p 1 l n - 1 ) h 0 l - p 0 l ( nh 1 l - 1 ) ( 1 - p 1 l ) h 0 l + ( h 1 l - 1 ) p 0 l + ( p 1 l - h 1 l ) , L 2 l = N 2 l = ( 1 - n ) h 0 l - 1 + nh 1 l ( 1 - p 1 l ) h 0 l + ( h 1 l - 1 ) p 0 l + ( p 1 l - h 1 l ) , L 3 l = N 3 l = ( n - 1 ) p 0 l + 1 - np 1 l ( 1 - p 1 l ) h 0 l + ( h 1 l - 1 ) p 0 l + ( p 1 l - h 1 l ) . The above solution is correct if and only if:
(72)
k + L 1 l + L 2 l p kl + L 3 l h kl = 0 , k = 2 , , n - 1 ;
(73)
( 1 - p 1 l ) h 0 l + ( - 1 + h 1 l ) p 0 l + p 1 l - h 1 l 0 , l = 1 , 2 .
Conditions (72)–(73) follow from (65), (66): both of the systems have 2 n + 2 equations and 6 unknowns. When n > 2 , the number of equations surpasses the number of unknown parameters, which means that the above conditions are necessary in order for the system to remain consistent.

As shown above, in the case of higher order solitary solutions, the auxiliary parameters L kl and N kl must be equal. Thus, inserting L kl = N kl ( k = 1 , 2 , 3 ) into (58) results in the following constraints on the Riccati system parameters:

(74)
a 2 = λ y ; b 2 = λ x ; c 2 = λ w ; d 2 = λ z .

Applying analogous reasoning used for the derivation of conditions (72)–(73) earlier, it could be concluded that solutions (70) to the system (67) hold if and only if the following relationships are true:

(75)
h k 2 + K 11 h k 1 + K 21 + K 31 kh k 1 = 0 , k = 1 , , n - 1 ,
(76)
K 1 l h k 2 + K 22 + K 32 kh k 2 = 0 , k = 1 , , n - 1 .

The results of the analysis performed in this subsection are summarized in the following Lemma.

Lemma 4.1

The Riccati system (2) admits higher order solitary solutions ( n 2 ) if and only if (72), (74) and (75)–(76) do hold true.

4.3

4.3 The algebraic structure of higher order solitary solutions ( n 2 )

The solutions (11)–(12) can be rewritten in a different form considering (61)–(62) and (64):

(77)
x t = x ^ t ^ = σ 1 X T = σ 1 t ^ n + p ( n - 1 ) 1 æ n - 1 t ^ n - 1 + + p 01 æ 0 t ^ n + æ n - 1 t ^ n - 1 + + æ 0 ;
(78)
y t = y ^ t ^ = σ 2 Y T = σ 2 t ^ n + h ( n - 1 ) 1 æ n - 1 t ^ n - 1 + + h 01 æ 0 t ^ n + æ n - 1 t ^ n - 1 + + æ 0 ,
(79)
z t = z ^ t ^ = σ 3 Z T = σ 3 t ^ n + p ( n - 1 ) 2 æ n - 1 t ^ n - 1 + + p 02 æ 0 t ^ n + æ n - 1 t ^ n - 1 + + æ 0 ;
(80)
w t = w ^ t ^ = σ 4 W T = σ 4 t ^ n + h ( n - 1 ) 2 æ n - 1 t ^ n - 1 + + h 02 æ 0 t ^ n + æ n - 1 t ^ n - 1 + + æ 0 .

It is important to note that necessary and sufficient conditions mentioned in Lemma 4.1 do result in a specific algebraic structure of higher order solitary solutions. The following Lemma splits the family of solutions into two cases.

Lemma 4.2

While the Riccati system (2) can admit solitary solutions of arbitrary order, some solutions may have constraints imposed on their parameters:

  • (i) The unconstrained case at n = 2 . The conditions (72), (75)–(76) are fulfilled without constraining the parameters of the solitary solution. Furthermore, any system can admit an infinite number of unique solitary solutions which can be generated by selecting different æ k , k = 0 , , n - 1 in (77)–(80).

  • (ii) The constrained case at n > 2 . The conditions (72), (75)–(76) are fulfilled only if n - 2 parameters æ k , k { 0 , , n - 1 } are set to zero. Selecting different values of the nonzero æ k , results in an infinite number of unique solitary solutions.

Note that case (ii) of the above Lemma does not inherently mean that higher order solitary solutions become more complex: due to the result that n - 2 parameters æ k , k { 0 , , n - 1 } must be set to zero, these solutions retain a maximum number of two extrema.

5

5 Numerical experiments

5.1

5.1 Kink solitary solution to a single Riccati equation

Consider the following Riccati equation:

(81)
x = 1 + 9 10 x - 5 8 x 2 ; x ( 0 ) = - 0.71 . The solution to (81) is a kink solitary wave depicted in Fig. 2. Regardless of initial conditions, the solutions to (81) are monotonous (Polyanin and Zaitsev, 2004).
Kink solitary solution x ( t ) to (81).
Fig. 2
Kink solitary solution x ( t ) to (81).

5.2

5.2 Bright/dark solitary solutions to a system of Riccati equations with multiplicative coupling

Consider the following system of Riccati equations with multiplicative coupling:

(82)
x = 1 + 9 10 x - 5 8 x 2 + 1 20 xy ;
(83)
y = - 33 10 + 3 10 y + 1 20 y 2 - 5 8 xy .
The solutions to (82)–(83) are bright/dark solitary waves depicted in Fig. 3. A detailed account of solitary solutions in such systems is given in Navickas et al. (2016). Note that compared to the single Riccati equation, the solitary solutions to this system are of a higher order and possess extrema.
Bright/dark solitary solutions x ( t ) (part (a)) and y ( t ) (part (b)) to (82)–(83).
Fig. 3
Bright/dark solitary solutions x ( t ) (part (a)) and y ( t ) (part (b)) to (82)–(83).

5.3

5.3 Second order solitary solutions to the system (2)

Suppose the following Riccati system is given:

(84)
x = 1 3 + 5869 3 x - 1730 3 x 2 - 460 xy ;
(85)
y = 275 18 + 105367 54 y - 460 y 2 - 1730 3 xy - 1 6 w - y ;
(86)
z = 1 + 9 10 z - 5 8 z 2 + 1 20 zw ;
(87)
w = - 33 10 + 3 10 w + 1 20 w 2 - 5 8 zw .
The techniques presented above are used to construct second order solitary solutions ( n = 2 ) (11)–(12).

Firstly, solitary solution parameters σ k ( k = 1 , 2 , 3 , 4 ), η , t 0 and æ k ( k = 0 , 1 ) are selected as follows:

(88)
σ 1 = 1 ; σ 2 = 3 ; σ 3 = 2 ; σ 4 = - 3 ; η = 1 ; t 0 = - 5 ; æ 0 = 1 ; æ 1 = 5 .

Then, the auxiliary parameters α l , β l , A kl , L ml , N ml and K ml ( m = 1 , 2 , 3 ; l , k = 1 , 2 ) are computed from (58):

(89)
α 1 = 5 3 ; α 2 = 3 2 ; β 1 = 5864 3 ; β 2 = - 3 5 ; A 11 = 5 ; A 12 = 3 ; A 21 = - 3 ; A 22 = - 2 ; L 11 = N 11 = 5864 3 ; L 12 = N 12 = - 3 5 ; L 21 = N 21 = - 1730 3 ; L 22 = N 22 = - 5 4 ; L 31 = N 31 = - 1380 ; L 32 = N 32 = - 3 20 ; K 11 = - 176 9 ; K 12 = - 27 10 ; K 21 = 275 9 ; K 22 = - 33 10 ; K 31 = - 6 ; K 32 = 3 .

Having derived the auxiliary parameters, it is now possible to use (61)–(62), (70), (71) and (75)–(76) to obtain the missing solitary solution parameters p ml and h ml ( m = 0 , 1 ; l = 1 , 2 ):

(90)
p 01 = - 1 5 ; p 02 = - 1 3 ; p 11 = - 1 2 ; p 12 = - 1 ; h 01 = 3 2 ; h 02 = - 11 9 ; h 11 = 187 115 ; h 12 = 11 .

The above parameters yield the following non-singular second order solitary solutions:

(91)
x t = 10 exp ( 2 t - 10 ) - 25 exp ( t - 5 ) - 2 10 exp ( 2 t - 10 ) + 5 exp ( t - 5 ) + 1 ;
(92)
y t = 3 46 exp ( 2 t - 10 ) + 69 exp ( t - 5 ) + 374 46 exp ( 2 t - 10 ) + 5 exp ( t - 5 ) + 1 ;
(93)
z t = 2 3 exp ( 2 t - 10 ) - 15 exp ( t - 5 ) - 1 3 exp ( 2 t - 10 ) + 5 exp ( t - 5 ) + 1 ;
(94)
w t = - 3 9 exp ( 2 t - 10 ) + 495 exp ( t - 5 ) - 11 9 exp ( 2 t - 10 ) + 5 exp ( t - 5 ) + 1 .
The above solutions are depicted in Fig. 4.
Second order solitary solutions x ( t ) (part (a)), y ( t ) (part (b)), z ( t ) (part (c)) and w ( t ) (part (d)) to (84)–(87).
Fig. 4
Second order solitary solutions x ( t ) (part (a)), y ( t ) (part (b)), z ( t ) (part (c)) and w ( t ) (part (d)) to (84)–(87).

5.4

5.4 Third order solitary solutions to the system (2)

In this subsection, the third order solitary solutions ( n = 3 ) are constructed for the following Riccati system:

(95)
x = 1 2 + 101 44 x - 30 11 x 2 - 1 44 xy ;
(96)
y = 15 4 - 5 11 y - 1 44 y 2 - 30 11 xy - w - y ;
(97)
z = 3 2 + z - 45 56 z 2 + 1 21 zw ;
(98)
w = - 15 2 - 3 4 w + 1 21 w 2 - 45 56 zw .

Analogously as in the previous subsection, solitary solution parameters σ k ( k = 1 , 2 , 3 , 4 ), η , t 0 and æ k ( k = 0 , 1 ) are selected as follows:

(99)
σ 1 = 1 ; σ 2 = 3 ; σ 3 = 2 ; σ 4 = - 3 ; η = 1 ; t 0 = - 5 ; æ 0 = 1 ; æ 1 = 10 .

From the Lemma 4.2 it follows that æ 2 = 0 .

The auxiliary parameters α l , β l , A kl , L ml , N ml and K ml ( m = 1 , 2 , 3 ; l , k = 1 , 2 ) are then obtained from (58):

(100)
α 1 = 5 2 ; α 2 = 9 4 ; β 1 = - 9 44 ; β 2 = - 5 4 ; A 11 = 5 ; A 12 = 3 ; A 21 = - 2 ; A 22 = - 4 3 ; L 11 = N 11 = - 9 44 ; L 12 = N 12 = - 5 4 ; L 21 = N 21 = - 30 11 ; L 22 = N 22 = - 45 28 ; L 31 = N 31 = - 3 44 ; L 32 = N 32 = - 1 7 ; K 11 = 3 4 ; K 12 = - 3 2 ; K 21 = 5 4 ; K 22 = - 15 2 ; K 31 = - 1 ; K 32 = 3 .

Finally, the solitary solution parameters p ml and h ml ( m = 0 , 1 ; l = 1 , 2 ) are computed using (61)–(62), (70), (71) and (75)–(76):

(101)
p 01 = - 1 5 ; p 02 = - 1 3 ; p 11 = - 1 3 ; p 12 = - 3 5 ; h 01 = 5 ; h 02 = - 5 ; h 11 = 25 ; h 12 = 5 .

Thus, the Riccati system (95)–(98) admits the following non-singular third order solitary solutions:

(102)
x t = 15 exp ( 3 t - 15 ) - 50 exp ( t - 5 ) - 3 15 exp ( 3 t - 15 ) + 10 exp ( t - 5 ) + 1 ;
(103)
y t = 3 exp ( 3 t - 10 ) + 250 exp ( t - 5 ) + 5 exp ( 3 t - 15 ) + 10 exp ( t - 5 ) + 1 ;
(104)
z t = 2 3 exp ( 3 t - 10 ) - 18 exp ( t - 5 ) - 1 3 exp ( 3 t - 15 ) + 10 exp ( t - 5 ) + 1 ;
(105)
w t = - 3 exp ( 3 t - 10 ) + 50 exp ( t - 5 ) - 5 exp ( 3 t - 15 ) + 10 exp ( t - 5 ) + 1 .
The solitary solutions (102)–(105) are displayed in Fig. 5. Note that although the order of the solitary solution is higher compared to the previous example (see Section 5.3), the number of extrema remains the same.
Third order solitary solutions x ( t ) (part (a)), y ( t ) (part (b)), z ( t ) (part (c)) and w ( t ) (part (d)) to (95)–(98).
Fig. 5
Third order solitary solutions x ( t ) (part (a)), y ( t ) (part (b)), z ( t ) (part (c)) and w ( t ) (part (d)) to (95)–(98).

6

6 Conclusions

Solitary solutions to model (2), representing cancer metastasis, were derived in this paper. Note that (2) represents a system with a one-way interaction between tumors: both tumor nodes are represented via Riccati systems with multiplicative coupling, however, the first node is autonomous (in the sense that it is not influenced by the remainder of the system), and the second is diffusively coupled with the first. This results in a system where the inter-node interactions are represented by a single diffusive term, rather than a two-way interaction consisting of two terms, discussed in detail in Telksnys et al. (2020).

While system (2) may appear similar to the Riccati system discussed in Telksnys et al. (2020), their dynamics are distinct. In either case, the evolution of tumor cell populations is non-monotonous, possessing both maxima and minima for finite time values as well as differing ratios of limits as time tends to - or + . However, with a single diffusive term, as discusssed in this paper, solitary solutions with at most two maxima (irregardles of the actual order n of the solution) can be obtained, while two diffusive terms in the analogous system results in at most three maxima.

Furthermore, it can be noted that the relation between cell population in the nodes themselves can be representative of tumor evolution: functions x ( t ) , z ( t ) evolve from a higher steady-state to a lower one, while y ( t ) , w ( t ) do the opposite. It is interesting to observe that the androgen-deprivation prostate cancer treatment model does admit only bright/dark (the second order) solitary solutions (Telksnys et al., 2020) However, this paper shows that the order of solitary solutions admitted by the metastasis model represented by multiplicatively and diffusively coupled Riccati systems can be larger than two. Together with the results of Telksnys et al. (2020), this paper provides a comprehensive study on solitary solutions to tumor models. Extensions of these models containing more complex coupling terms (both inside the nodes and between them) remain a definite objective of future research.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

References

  1. Alqudah, M.A., 2020. Cancer treatment by stem cells and chemotherapy as a mathematical model with numerical simulations, Alexandria Eng. J.
  2. , , , . The mathematics of cancer: integrating quantitative models. Nat. Rev. Cancer. 2015;15(12):730-745.
    [Google Scholar]
  3. , , , , . Maximization of viability time in a mathematical model of cancer therapy. Mathe. Biosci.. 2017;294:110-119.
    [Google Scholar]
  4. , , , , , . Multiple primary cancers: An enigma. South Asian J. Cancer. 2016;5(1):29.
    [Google Scholar]
  5. , , , , , , , , , , . Clinicopathologic features of metachronous or synchronous gastric cancer patients with three or more primary sites. Cancer Res. Treat.: Off. J. Korean Cancer Assoc.. 2010;42(4):217.
    [Google Scholar]
  6. , , , , , , , , . Prediction of metachronous multiple primary cancers following the curative resection of gastric cancer. BMC Cancer. 2013;13(1):394.
    [Google Scholar]
  7. , , . A mathematical model of cancer treatment by radiotherapy followed by chemotherapy. Mathe. Comput. Simul.. 2016;124:1-15.
    [Google Scholar]
  8. , , , . Current status of mathematical modeling of cancer–from the viewpoint of cancer hallmarks. Curr. Opin. Syst. Biol.. 2017;2:39-48.
    [Google Scholar]
  9. , . Mathematical modeling of cancer metabolism. Crit. Rev. Oncol./Hematol.. 2018;124:37-40.
    [Google Scholar]
  10. , . Multiple primary malignant neoplasms. historical perspectives. Cancer. 1977;40(S4):1786-1792.
    [Google Scholar]
  11. , . Migration difference in diffusively-coupled prey–predator system on heterogeneous graphs. Physica A. 2020;537:122705
    [Google Scholar]
  12. , , , . Existence of solitary solutions in a class of nonlinear differential equations with polynomial nonlinearity. Appl. Math. Comput.. 2016;283:333-338.
    [Google Scholar]
  13. , , , , . Existence of second order solitary solutions to riccati differential equations coupled with a multiplicative term. IMA J. Appl. Mathe.. 2016;81(6):1163-1190.
    [Google Scholar]
  14. , , . Handbook of Nonlinear Partial Differential Equations. CRC Press; .
  15. , , , . Enhancing master-slave synchronization: The effect of using a dynamic coupling. Phys. Rev. E. 2018;98(1):012208
    [Google Scholar]
  16. , , . A mathematical model for the immune-mediated theory of metastasis. J. Theoret. Biol.. 2019;482:109999
    [Google Scholar]
  17. , , . Application of top-down and bottom-up systems approaches in ruminant physiology and metabolism. Curr. Genomics. 2012;13(5):379-394.
    [Google Scholar]
  18. , , , , . Mathematical modelling of cancer stem cell-targeted immunotherapy. Math. Biosci.. 2019;318:108269
    [Google Scholar]
  19. , , , , , . Environmental factors in breast cancer invasion: a mathematical modelling review. Pathology. 2017;49(2):172-180.
    [Google Scholar]
  20. , , , , , . Homoclinic and heteroclinic solutions to a hepatitis c evolution model. Open Mathe.. 2018;16(1):1537-1555.
    [Google Scholar]
  21. , , , , , . Symmetry breaking in solitary solutions to the hodgkin–huxley model. Nonlinear Dyn.. 2019;97(1):571-582.
    [Google Scholar]
  22. , , , , , , . Solitary solutions to an androgen-deprivation prostate cancer treatment model. Mathe. Methods Appl. Sci.. 2020;43(7):3995-4006.
    [Google Scholar]
  23. , , . Analysis of mathematical model of prostate cancer with androgen deprivation therapy. Commun. Nonlinear Sci. Numer. Simul.. 2019;66:41-60.
    [Google Scholar]
Show Sections