Next Article in Journal
Lateral Fractal Formation by Crystallographic Silicon Micromachining
Next Article in Special Issue
Dynamic Properties for a Second-Order Stochastic SEIR Model with Infectivity in Incubation Period and Homestead-Isolation of the Susceptible Population
Previous Article in Journal
Differential Subordination and Superordination Results for q-Analogue of Multiplier Transformation
Previous Article in Special Issue
A Compact Scheme Combining the Fast Time Stepping Method for Solving 2D Fractional Subdiffusion Equations
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Dynamics and Stability of a Fractional-Order Tumor–Immune Interaction Model with B-D Functional Response and Immunotherapy

1
School of Science, Xi’an Technological University, Xi’an 710032, China
2
School of Mathematics and Statistics, Xi’an Jiaotong University, Xi’an 710049, China
*
Author to whom correspondence should be addressed.
Fractal Fract. 2023, 7(2), 200; https://doi.org/10.3390/fractalfract7020200
Submission received: 31 December 2022 / Revised: 30 January 2023 / Accepted: 4 February 2023 / Published: 17 February 2023
(This article belongs to the Special Issue Fractional Differential Equations in Anomalous Diffusion)

Abstract

:
In this paper, we investigate a fractional-order tumor–immune interaction model with B-D function item and immunotherapy. First, the existence, uniqueness and nonnegativity of the solutions of the model are established. Second, the local and global asymptotic stability of some tumor-free equilibrium points and a unique positive equilibrium point are obtained. Finally, we use numerical simulation method to visualize and verify the theoretical conclusions. It is known that the fractional-order parameter β has a stabilization effect, and the tumor cells can be destroyed or controlled by using immunotherapy.

1. Introduction

In this paper, we mainly investigated a fractional-order tumor–immune interaction model with B-D function item and immunotherapy as follows
{ D 0 c t β u ( t ) = a v m 1 u + k 1 u w 1 + w + s 1 , D 0 c t β v ( t ) = v ( b v ) c u v u + v + d , D 0 c t β w ( t ) = k 2 u v u + v + d m 2 w + s 2 , u ( 0 ) = u 0 0 , v ( 0 ) = v 0 0 , w ( 0 ) = w 0 0 ,
where u ( t ) stand for the density of effector cells that are cytotoxic to the tumor cells; v ( t ) represents the density of tumor cells; w ( t ) is the concentration of interleukin-2 (IL-2); a is the antigenicity of the tumor; b is the environmental carrying capacity; k 1 ,   c ,   k 2 are the uptake velocity when all sites are saturated by the substrate; 1 / m 1 is the natural average lifespan; m 2 is the loss or degraded rate of IL-2; d is the half saturation constant; s 1 represents the treatment by an external source of effector cells (ACI); s 2 represents the treatment by an external input of IL-2 into the model, D 0 c t β is the standard Caputo differentiation; and u v / ( u + v + d ) is the Beddington–DeAngelis (denoted by B-D) functional response.
A tumor is a new organism formed by local tissue cell proliferation under the action of various tumorigenic factors. Tumors can also be classified as benign, precancerous, or malignant, with malignant tumors also known as cancers. Because cancer cells divide uncontrollably and can invade other tissues, the individual fatality rate for cancer is extremely high. The safe treatment of cancer is an important issue being studied in the medical model at present. Numerous studies have shown that host immune cells have an inhibitory effect on tumor growth [1,2,3,4,5,6,7].
Immunotherapy is a cancer treatment that utilizes cytokines and adoptive cell immunotherapy (ACI). Moreover, interleukin-2 (IL-2), the major cytokine responsible for lymphocyte activation, growth, and differentiation, is mainly produced by CD4+ T cells. The injection of cultured immune cells with antitumor reactivity into a tumor-carrying host is performed by ACI, using either lymphokine-activated killer (LAK) therapy or tumor-infiltrating lymphocyte (TIL) therapy in combination with abundant IL-2. By applying LAK therapy or TIL therapy or by applying both therapies simultaneously, Krischner et al. [8] researched a model describing tumor–immune dynamics together with the feature of IL-2 dynamics. In [9], we have proposed a fractional-order model with a Holling II function item describing the interaction between effector cells, tumor cells, and IL-2 as follows
{ D 0 c t β u ( t ) = a v m 1 u + k 1 u w 1 + w + s 1 , D 0 c t β v ( t ) = v ( 1 b v ) c u v v + d , D 0 c t β w ( t ) = k 2 u v 1 + v m 2 w + s 2 , u ( 0 ) = u 0 0 , v ( 0 ) = v 0 0 , w ( 0 ) = w 0 0 ,
where D 0 c t β is the standard Caputo differentiation and β ( 0 ,   1 ) , the Caputo fractional derivative of order β , is defined as [10,11]
D 0 c t β h ( x ) = 1 Γ ( m β ) t 0 t ( x t ) m β 1 h ( m ) ( t ) d t ,   m 1 < β < m ,   m .
Fractional-order differential equations present an outstanding description of certain nonlinear practical problems. Many real biological models have memory, and fractional models may be better suited to describing the real situation than integer models, which have the advantage that they describe the entire time domain of a physical process, whereas integer models relate to the local properties of a location. So, we have studied the asymptotic behavior analysis of a fractional-order tumor–immune interaction model with Holling functional response function and immunotherapy in [9].
However, functional responses are of important biological significance in the model (Equation (2)), which represents the average conversion rate among three types of cells. It is easy to see that the above Holling II reaction functions in (Equation (2)) are only dependent on the transformed cells and are not affected by the number of the newborn cells. In order to formulate mutual interference among cells, Beddington [12] and DeAngelis [13] established so-called Beddington–DeAngelis (denoted by B-D) functional response given by u v u + v + d . Therefore, by introducing the B-D functional response into the model (Equation (2)) and taking into account the fractional-order equation and the important influence of effector cells on the immunotherapy process, we established model (Equation (1)), which can more accurately simulate the transformation relations and laws between effector cells, tumor cells, and IL-2, and enrich the research results of fractional biological mathematical ODE models.
This paper is organized as follows. In Section 2, we state our main results about the existence, uniqueness, and nonnegativity of the solutions of the model (Equation (1)) by using the standard comparison theorem and the local Lipschitz condition. Meanwhile, there are some equilibrium points, depending on what treatment is not used and whether the tumor is curable. Additionally, sufficient conditions on the local and global asymptotic stability of these equilibrium points are proved using the linearization method and the Routh–Hurwitz conditions in Section 3, for the cases of the realistic tumor-free equilibriums and the unique interior equilibriums, respectively. Finally, some numerical simulations are given to illustrate and visualize our theoretical results in Section 4.

2. Uniqueness and Nonnegativity of Solutions

In this section, we study the existence, uniqueness, and nonnegativity of the solutions of the model (Equation (1)). Then, we need the following lemma to prove the existence and uniqueness of the solutions of the fractional-order model (Equation (1)).
Lemma 1
(see [14,15]). Consider the model with initial condition  h ( t 0 )
D t 0 c t β h ( t ) = h ( t , x ) ,
where β ( 0 , 1 ) and  h : [ t 0 , ) × Ω N , Ω N , if  h ( t , x )  meets the local Lipschitz condition depending on  x , then the model presents a unique solution of Equation (3).
Theorem 1.
Define  Ω = { ( u , v , w ) 3 : max { | u | , | v | , | w | } M ˜ } . For every initial condition X 0 =   ( u 0 , v 0 , w 0 ) Ω , the model (Equation (1)) has a unique solution for all  t 0 .
Proof. 
Set T ( 0 , + ) ; we denote X = ( u , v , w ) T and X ˜ = ( u ˜ , v ˜ , w ˜ ) T in the region Ω × ( 0 , T ] . Consider a mapping defined by H ( X ) = ( H 1 ( X ) , H 2 ( X ) , H 3 ( X ) ) T , where
H 1 ( X ) = a v m 1 u + k 1 u w 1 + w + s 1 , H 2 ( X ) = v ( b v ) c u v u + v + d , H 3 ( X ) = k 2 u v u + v + d m 2 w + s 2 .
For any X , X ˜ Ω , according to (Equation (4)), we obtain the following inequality
H ( X ) H ( X ˜ ) = | H 1 ( X ) H 1 ( X ˜ ) | + | H 2 ( X ) H 2 ( X ˜ ) | + | H 3 ( X ) H 3 ( X ˜ ) | = | a v m 1 u + k 1 u w 1 + w + s 1 a v ˜ + m 1 u ˜ k 1 u ˜ w ˜ 1 + w ˜ s 1 | + | v ( b v ) c u v u + v + d v ˜ ( b v ˜ ) + c u ˜ v ˜ u ˜ + v ˜ + d | + | k 2 u v u + v + d m 2 w + s 2 k 2 u ˜ v ˜ u ˜ + v ˜ + d + m 2 w ˜ s 2 | = | a ( v v ˜ ) + ( k 1 w ˜ 1 + w ˜ m 1 ) ( u u ˜ ) + k 1 u ˜ ( w w ˜ ) ( 1 + w ) ( 1 + w ˜ ) | + | [ b ( v + v ˜ ) c u u + v + d + c u v ˜ ( u + v + d ) ( u ˜ + v ˜ + d ) ] ( v v ˜ ) c v ˜ ( v + d ) ( u + v + d ) ( u ˜ + v ˜ + d ) ( u u ˜ ) | + | k 2 u ˜ ( u + d ) ( u + v + d ) ( u ˜ + v ˜ + d ) ( v v ˜ ) + [ k 2 v u + v + d k 2 u ˜ v ( u + v + d ) ( u ˜ + v ˜ + d ) ] ( u u ˜ ) m 2 ( w w ˜ ) | a | v v ˜ | + ( k 1 + m 1 ) | u u ˜ | + k 1 M ˜ | w w ˜ | + 2 k 2 | u u ˜ | + k 2 | v v ˜ | + m 2 | w w ˜ | +   ( b + 2 c + 2 M ˜ ) | v v ˜ | + c | u u ˜ | = ( k 1 + m 1 + 2 k 2 + c ) | u u ˜ | + ( a + k 2 + b + 2 c + 2 M ˜ ) | v v ˜ | + ( k 1 M ˜ + m 2 ) | w w ˜ | L X X ˜ ,
where L = max { k 1 + m 1 + 2 k 2 + c , a + k 2 + b + 2 c + 2 M ˜ , k 1 M ˜ + m 2 } ; then H ( X ) with respect to X satisfies the local Lipschitz condition, and model (Equation (3)) has a unique solution.
Thanks to the above analysis, in the following, we establish the conditions for the nonnegativity of the solutions of the model (Equation (1)). □
Theorem 2.
We denote  * = { x : x 0 }  and  Ω * = { ( u , v , w ) Ω : u , v   a n d   w * } . For every initial condition  X 0 = ( u 0 , v 0 , w 0 ) + 3 , all solutions of model (Equation (1)) are nonnegative.
Proof. 
To end this Theorem, we will adopt the method of contradiction proof. Assume that there exists t ˜ [ 0 , + ) such that the solutions of model (Equation (1)) pass through the u , v , or w axes. Set β to be defined in (Equation (3)); the results can be divided into the following three cases:
(i)
If u ( t ˜ ) > 0 , v ( t ˜ ) > 0 , and w ( t ˜ ) = 0 , there is t ˜ > t ˜ and 0 < t ˜ t ˜ 1 . When t [ t ˜ , t ˜ ) , it is known that u ( t ) > 0 , v ( t ) > 0 and w ( t ) < 0 . From the third equation of model (Equation (1)), we can derive D 0 c t β w ( t ) > m 2 w and then w ( t ) > w ( 0 ) E β [ ( m 2 ) t β ] with t [ t ˜ , t ˜ ) . The Mittag–Leffler function is E β , γ ( t ) = s = 0 [ t s / Γ ( s β + γ ) ] and E β ( t ) = E β , 1 ( t ) with β , γ > 0 [16]. We know that the Mittag–Leffler function is positive in [11,14,15], then E β ( t ) > 0 . By using the standard comparison theorem [9,17], we obtain w ( t ) > 0 for all t [ t ˜ , t ˜ ) , and β ( 0 , 1 ) , which contradicts the hypothesis. The same can be said for the following two cases:
(ii)
Suppose that u ( t ˜ ) > 0 , w ( t ˜ ) > 0 and v ( t ˜ ) = 0 , there is t ˜ > t ˜ and 0 < t ˜ t ˜ 1 . When t [ t ˜ , t ˜ ) , it is known that u ( t ) > 0 , w ( t ) > 0 , and v ( t ) < 0 . From the second equation of the model (Equation (1)), we can derive D 0 c t β v ( t ) > b v and then v ( t ) > v ( 0 ) E β ( b t β ) with t [ t ˜ , t ˜ ) . Then, we obtain v ( t ) > 0 for all t [ t ˜ , t ˜ ) , and β ( 0 , 1 ) , which contradicts the hypothesis.
(iii)
Suppose that v ( t ˜ ) > 0 , w ( t ˜ ) > 0 , and u ( t ˜ ) = 0 , there is t ˜ > t ˜ and 0 < t ˜ t ˜ 1 . When t [ t ˜ , t ˜ ) , it is known that v ( t ) > 0 , w ( t ) > 0 , and u ( t ) < 0 . From the first equation of model (Equation (1)), we can derive D 0 c t β u ( t ) > ( k 1 m 1 ) u and then u ( t ) > u ( 0 ) E β [ ( k 1 m 1 ) t β ] with t [ t ˜ , t ˜ ) . Then, we obtain u ( t ) > 0 for all t [ t ˜ , t ˜ ) , and β ( 0 , 1 ) , which contradicts the hypothesis. □
According to the above method, the nonnegativity of the other solutions can be proved, and the proof process is omitted here.
According to the above cases, it follows that the solutions of model (Equation (1)) are nonnegative.

3. Equilibria and Stability

In this section, all positive constant equilibrium points of model (Equation (1)) will be analyzed. From a direct calculation, we know that model (Equation (1)) has four nonnegative equilibrium points in order to simulate the disappearance of the tumor cells; all the points at least have one zero element.
(i)
A 1 = ( 0 , 0 , 0 ) , if s 1 = 0 and s 2 = 0 ,
(ii)
A 2 = ( s 1 m 1 , 0 , 0 ) , if s 1 > 0 and s 2 = 0 , tumor treatment by an external source of effector cells,
(iii)
A 3 = ( 0 , 0 , s 2 m 2 ) , if s 1 = 0 and s 2 > 0 ,
(iv)
A 4 = ( s 1 ( m 2 + s 2 ) m 1 m 2 + s 2 ( m 1 k 1 ) ,   0 ,   s 2 m 2 ) , if s 1 > 0 , s 2 > 0 and m 1 m 2 > s 2 ( k 1 m 1 ) , tumor treatment by an external source of effector cells with an external input of IL-2 into the model (Equation (1)).
The points (i) and (iii) are not in line with the scope of actual medical research, because effector cells do not disappear under any circumstances; thus, these two points will not be discussed in depth in this section. However, the cases (ii) and (iv) are tumor-free equilibrium points in reality, it remains that if the tumor is cured, we will study these two cases and provide sufficient conditions, which is necessary for asymptotic stability in cases (ii) and (iv). Then, we only give the conditions that make model (Equation (1)) describe a unique positive equilibrium point A ^ = ( u ^ , v ^ , w ^ ) and omit calculation the procedure. We can easily obtain the following results from [5].
Lemma 2.
If one of the following inequalities is true, model (Equation (1)) has a unique positive equilibrium point  A ^ .
s 1 > a ( d c ) ,   ( c b d ) s 1 2 + ( a b c a c d a b 2 + a d 2 ) s 1 < a 2 b d ( c d ) ,   ( m 1 k 1 ) b d c b + v > s 1   and   m 1 k 1 .
s 1 > a d ,   c s 1 2 + ( a b c + a d 2 ) s 1 < a 2 b d ( c d ) ,   m 1 b d c b > s 1 ( 1 + s 2 )   and   m 1 = k 1 .
where A ^ = ( u ^ , v ^ , w ^ ) is given as follows
u ^ = v ^ 2 + ( b s ) v ^ + b d c b + v ^ , w ^ = m 1 u ^ a v ^ s 1 ( k 1 m 1 ) u ^ + a v ^ + s 1 , w ^ = ( s 2 + k 2 v ^ ) u ^ + ( v ^ + d ) s 2 u ^ + v ^ + d .
In the following lemma, we give the conditions for the asymptotic stability.
Lemma 3.
Let  J ( X )  be the Jacobian matrix of model (Equation (1)); we linearize model (Equation (1)) at point  X = ( u , v , w )  and obtain
J ( X ) = H X ( X ) = ( m 1 + k 1 w 1 + w a k 1 u ( 1 + w ) 2 c v ( v + d ) ( u + v + d ) 2 b 2 v c u ( u + d ) ( u + v + d ) 2 0 k 2 v ( v + d ) ( u + v + d ) 2 k 2 u ( u + d ) ( u + v + d ) 2 m 2 ) ,
where  H ( X )   is defined in Theorem 1, let   λ j ( j = 1 , 2 , ) be the eigenvalues of J ( X ) . If all eigenvalues λ j ( j = 1 , 2 , ) satisfy | arg ( λ i ) | > β π 2 , equilibrium point X is a local asymptotic stability; if some eigenvalues λ j ( j = 1 , 2 , ) satisfy | arg ( λ i ) | > β π 2 and others satisfy | arg ( λ i ) | < β π 2 , X   is a saddle point and unstable.
In the following, we give some conditions to prove the local asymptotic stability of A 2 and A 4 defined by Lemma 3.
Theorem 3.
Let  b ,   c , d , m 1 , and s 1  be defined in model (Equation (1)); we obtain the following statements
(i) 
If ( c b ) s 1 > m 1 b d , then equilibrium point A 2 is a local asymptotic stability.
(ii) 
If ( c b ) s 1 < m 1 b d , then equilibrium point A 2 is a saddle point and unstable.
Proof. 
Substituting X = A 2 into (Equation (7)), we have
J ( A 2 ) = ( m 1 a k 1 s 1 m 1 0 b c s 1 s 1 + m 1 d 0 0 k 2 s 1 s 1 + m 1 d m 2 ) ,
it is easy to see that the eigenvalues of J ( A 2 ) are
λ 1 = m 1 , λ 2 = m 2 , λ 3 = b c s 1 s 1 + m 1 d ,
if ( c b ) s 1 > m 1 b d , that is arg ( λ 3 ) = π , then all the other λ i ( i = 1 , 2 , 3 ) will also satisfy arg ( λ 1 ) = arg ( λ 2 ) = arg ( λ 3 ) = π . Therefore, | arg ( λ i ) | > β π 2 ( i = 1 , 2 , 3 ) for β ( 0 , 1 ) . It follows from Lemma 2 that the equilibrium point A 2 is a local asymptotic stability. However, if ( c b ) s 1 < m 1 b d , arg ( λ 1 ) = arg ( λ 2 ) = π and arg ( λ 3 ) = 0 , based on Lemma 2, | arg ( λ i ) | < β π 2 ( i = 1 , 2 , 3 ) for β ( 0 , 1 ) . Then, the equilibrium point A 2 is a saddle point and unstable. □
Theorem 4.
Let  b ,   c , d , m 1 ,   m 2 ,   s 1 , and s 2  be defined in model (Equation (1)); we obtain the following statements
(i)
If s 2 < m 1 ( m 2 + s 2 ) and m 1 m 2 b d + s 2 b d ( m 1 k 1 ) < ( c s 1 b s 1 ) ( m 2 + s 2 ) , then equilibrium A 4 is a local asymptotic stability.
(ii)
If s 2 > m 1 ( m 2 + s 2 ) or m 1 m 2 b d + s 2 b d ( m 1 k 1 ) > ( c s 1 b s 1 ) ( m 2 + s 2 ) , then equilibrium A 2 is a saddle point and unstable.
Proof. 
Substituting X = A 4 into (Equation (7)), we obtain
J ( A 4 ) = ( m 1 + s 2 m 2 + s 2 a ϒ 0 b c s 1 ( m 2 + s 2 ) s 1 ( m 2 + s 2 ) + m 1 m 2 d + s 2 d ( m 1 k 1 ) 0 0 k 2 s 1 ( m 2 + s 2 ) s 1 ( m 2 + s 2 ) + m 1 m 2 d + s 2 d ( m 1 k 1 ) m 2 ) ,
where ϒ = m 2 m 2 + s 2 k 1 s 1 m 1 m 2 + ( m 1 k 1 ) s 2 and it is clear that the eigenvalues of J ( A 4 ) are
λ 1 = m 2 ,
λ 2 = m 1 + s 2 m 2 + s 2 ,
λ 3 = b c s 1 ( m 2 + s 2 ) s 1 ( m 2 + s 2 ) + m 1 m 2 d + s 2 d ( m 1 k 1 ) .
If s 2 < m 1 ( m 2 + s 2 ) and m 1 m 2 b d + s 2 b d ( m 1 k 1 ) < ( c s 1 b s 1 ) ( m 2 + s 2 ) , that is arg ( λ 2 ) = arg ( λ 3 ) = π , then all the other λ i ( i = 1 , 2 , 3 ) will also satisfy arg ( λ 1 ) = arg ( λ 2 ) = arg ( λ 3 ) = π . Therefore, from Lemma 2, we obtain | arg ( λ i ) | >   β π 2 ( i = 1 , 2 , 3 ) for β ( 0 , 1 ) , the equilibrium A 4 is a local asymptotic stability. However, if s 2 > m 1 ( m 2 + s 2 ) , arg ( λ 2 ) = 0 . Additionally, if m 1 m 2 b d + s 2 b d ( m 1 k 1 ) > ( c s 1 b s 1 ) ( m 2 + s 2 ) , arg ( λ 3 ) = 0 , then the equilibrium point A 4 is a saddle point and unstable. □
In the following, applying the Hurwitz conditions, we roughly provide the local asymptotically stability of the unique positive equilibrium point A ^ .
Theorem 5.
Let  B ( q ) and q 3  be defined in (Equations (8) and (9)); if  m 1 k 1 b + m 2 > 0 , q 3 > 0 and B ( q ) > 0 , then A ^  is locally asymptotically stable.
Proof. 
Substituting X = A ^ into (Equation (7)), we have
J ( A ^ ) = ( m 1 + k 1 w ^ 1 + w ^ a k 1 u ^ ( 1 + w ^ ) 2 c v ^ ( v ^ + d ) ( u ^ + v ^ + d ) 2 b 2 v ^ c u ^ ( u ^ + d ) ( u ^ + v ^ + d ) 2 0 k 2 v ^ ( v ^ + d ) ( u ^ + v ^ + d ) 2 k 2 u ^ ( u ^ + d ) ( u ^ + v ^ + d ) 2 m 2 ) ,
treating J ( A ^ ) as the roots of K ( λ ) = λ 3 + q 1 λ 2 + q 2 λ + q 3 = 0 , then we can obtain the following equations
q 1 = m 1 k 1 w ^ 1 + w ^ b + 2 v ^ + c u ^ ( u ^ + d ) ( u ^ + v ^ + d ) 2 + m 2 , q 2 = ( k 1 w ^ 1 + w ^ m 1 ) [ b 2 v ^ m 2 c u ^ ( u ^ + d ) ( u ^ + v ^ + d ) 2 ] m 2 [ b 2 v ^ c u ^ ( u ^ + d ) ( u ^ + v ^ + d ) 2 ] + a c v ^ ( v ^ + d ) ( u ^ + v ^ + d ) 2 k 1 u ^ ( 1 + w ^ ) 2 k 2 v ^ ( v ^ + d ) ( u ^ + v ^ + d ) 2 , q 3 = [ b 2 v ^ c u ^ ( u ^ + d ) ( u ^ + v ^ + d ) 2 ] [ k 1 u ^ ( 1 + w ^ ) 2 k 2 v ^ ( v ^ + d ) ( u ^ + v ^ + d ) 2 m 1 m 2 + k 1 m 2 w ^ 1 + w ^ ] + c v ^ ( v ^ + d ) ( u ^ + v ^ + d ) 2 k 2 u ^ ( u ^ + d ) ( u ^ + v ^ + d ) 2 k 1 u ^ ( 1 + w ^ ) 2 + a c m 2 v ^ ( v ^ + d ) ( u ^ + v ^ + d ) 2 ,
if
q 1 = m 1 k 1 w ^ 1 + w ^ b + 2 v ^ + c u ^ ( u ^ + d ) ( u ^ + v ^ + d ) 2 + m 2 > m 1 k 1 b + m 2 > 0 , q 3 > 0 ,
B ( q ) = q 1 q 2 q 3 > 0
according to the Hurwitz conditions [18], all eigenvalues of J ( A ^ ) have negative real parts, we obtain the local stability conditions of A ^ . □
Next, we present the conditions for the global stability of A ^ .
Theorem 6.
Suppose that  m 1 > k 1 . Then, the nonnegative solution ( u , v , w ) to (Equation (1)) satisfies
lim t sup   Ω u ( x , t ) a b + s 1 m 1 k 1 , lim t sup   Ω v ( x , t ) b , lim t sup   Ω w ( x , t ) 1 m 2 ( k 2 a b + s 1 m 1 k 1 + s 2 ) .
Proof. 
To end this Theorem, we will adopt the method of the comparison principle [17]. From the second equation of the model (Equation (1)), we can verify that the sufficient conditions v ( b v ) c u v u + v + d v ( b v ) . Then, there exists some t 1 ( 0 , ) such that v ( x , t ) b + ε 1 for any positive constant ε 1 . Notice with this result, we obtain that a v m 1 u + k 1 u w 1 + w + s 1 s 1 ( m 1 k 1 ) u + a b with [ t 1 , ) × Ω ; thus, there exists t 2 [ t 1 , ) such that u ( x , t ) a b + s 1 m 1 k 1 + ε 2 for any nonnegative ε 2 . Finally, we obtain directly that there exists some t 3 [ t 2 , ) such that β 1 m 2 ( k 2 a b + s 1 m 1 k 1 + s 2 ) + ε 3 for any constant ε 3 on [ t 2 , ) × Ω . Therefore, by the arbitrariness of ε 1 , ε 2 , and ε 3 , the desired result has been obtained. □
Next, we use the Lyapunov function method to study the global stability of the unique positive equilibrium point A ^ . We first define the following Lyapunov function
S ( t ) S ( u ( t ) , v ( t ) , w ( t ) ) = 1 2 ( u u ^ ) 2 + ( v v ^ v ^ ln v v ^ ) + 1 2 ( w w ^ ) 2 ,
thanks to the famous Lyapunov stability; if D 0 c t β S ( t ) 0 , the unique positive equilibrium A ^ is global stability.
Lemma 4
(see [19]). Set  x ( t ) * and β ( 0 , 1 ) , we obtain the following inequality and  t > t 0 .
D 0 c t β [ x ( t ) x ^ x ^ ln x ( t ) x ^ ] [ 1 x ^ x ( t ) ] D 0 c t β x ( t ) ,   x ^ * .
Lemma 5
(see [20]). Let  x ( t ) *  be defined as a continuous and derivable function. Then, we obtain the following inequality and  t > t 0
1 2 D 0 c t β x 2 ( t ) x ( t ) D 0 c t β x ( t ) .
Theorem 7.
Assume that (Equation (5)) holds. Then, the positive constant solution  A ^  to (Equation (1)) is globally asymptotically stable if
{ m 1 + k 1 + a 2 + c d + k 2 + k 1 2 a b + s 1 m 1 k 1 < 0 , 2 c d 1 + a 2 + k 2 2 < 0 , m 2 + 3 2 k 2 + k 1 2 a b + s 1 m 1 k 1 < 0
.
Proof. 
To end this Theorem, we will use the model (Equation (1)) and integrate by parts, we have
D 0 c t β S ( t ) ( u u ^ ) D 0 c t β u ( t ) + [ 1 v ^ v ( t ) ] D 0 c t β v ( t ) + ( w w ^ ) D 0 c t β w ( t ) = ( u u ^ ) ( a v m 1 u + k 1 u w 1 + w + s 1 ) + ( v v ^ ) ( b v c u u + v + d ) + ( w w ^ ) ( k 2 u v u + v + d m 2 w + s 2 ) .
According to the definition of A ^ , Lemma 4 and Lemma 5, we derive
D 0 c t β S ( t ) ( u u ^ ) D 0 c t β u ( t ) + [ 1 v ^ v ( t ) ] D 0 c t β v ( t ) + ( w w ^ ) D 0 c t β w ( t ) = ( u u ^ ) ( a v m 1 u + k 1 u w 1 + w + s 1 ) + ( v v ^ ) ( b v c u u + v + d ) + ( w w ^ ) ( k 2 u v u + v + d m 2 w + s 2 ) = ( u u ^ ) 2 ( m 1 + k 1 w 1 + w ) + ( v v ^ ) 2 ( c u ^ ( u + v + d ) ( u ^ + v ^ + d ) 1 ) + ( u u ^ ) ( v v ^ ) ( a + c u ^ ( u + v + d ) ( u ^ + v ^ + d ) c u + v + d ) + k 2 u ^ ( u + d ) ( v v ^ ) ( w w ^ ) ( u + v + d ) ( u ^ + v ^ + d ) m 2 ( w w ^ ) 2 + ( u u ^ ) ( w w ^ ) ( k 1 u ^ ( 1 + w ) ( 1 + w ^ ) + k 2 v u + v + d k 2 u ^ v ( u + v + d ) ( u ^ + v ^ + d ) ) H ( t ) ,
Thanks to Theorem 5, it is easy to see that u ( x , t ) a b + s 1 m 1 k 1 , then
H ( t ) ( u u ^ ) 2 ( m 1 + k 1 ) + ( u u ^ ) ( v v ^ ) ( a + 2 c d ) + ( v v ^ ) 2 ( c d 1 ) + k 2 ( v v ^ ) ( w w ^ ) m 2 ( w w ^ ) 2 + ( u u ^ ) ( w w ^ ) ( k 1 u ^ + 2 k 2 ) ( u u ^ ) 2 ( m 1 + k 1 + a 2 + c d + k 2 + k 1 2 a b + s 1 m 1 k 1 ) + ( v v ^ ) 2 ( c d 1 + a 2 + c d + k 2 2 ) + ( w w ^ ) 2 ( m 2 + 3 2 k 2 + k 1 2 a b + s 1 m 1 k 1 ) .
By direct scaling calculation, we easily prove that D 0 c t β S ( t ) 0 and the positive constant solution A ^ to (Equation (1)) is globally asymptotically stable if the following three conditions are satisfied
(i)
m 1 + k 1 + a 2 + c d + k 2 + k 1 2 a b + s 1 m 1 k 1 < 0 ,
(ii)
2 c d 1 + a 2 + k 2 2 < 0 ,
(iii)
m 2 + 3 2 k 2 + k 1 2 a b + s 1 m 1 k 1 < 0 .
So, we complete the proof of Theorem 7. □

4. Numerical Simulations

In this section, we will obtain some numerical simulations by using the fractional linear multistep methods and the predictor–corrector algorithm [21,22], which are applied to find the approximate solution of the fractional-order tumor–immune interaction model (Equation (1)). Therefore, we can visualize and verify the above theoretical conclusions. Note that from the above numerical simulation results, we easily obtain the following conclusions.
(i)
From Theorem 3, it follows that when s 1 = 3 and s 2 = 0 , there is a realistic tumor-free equilibrium A 2 = ( 3 , 0 , 0 ) with a = 0.9 , m 1 = 1 , k 1 = 0.5 , b = 1 , c = 3 , s 1 = 3 , d = 2.5 , k 2 = 1 , m 2 = 1 , s 2 = 0 , and β = 0.9 . It can be observed in Figure 1 that the tumor-free equilibrium A 2 = ( 3 , 0 , 0 ) is locally asymptotically stable and it remains that the tumor is curable by using ACI therapy or ACI therapy and IL-2.
(ii)
According to Theorem 4, we take s 1 = 3 and s 2 = 1 , and there exists the tumor-free equilibrium A 4 = ( 4 , 0 , 1 ) . See Figure 2; A 4 is locally asymptotically stable and remains that the tumor is curable by ACI therapy and the external input of IL-2 into the model (Equation (1)).
(iii)
From Figure 3, it is easy to see that as tumor cells v of A 2 and A 4 are destroyed more quickly, the greater value of β . It is clearly known that the fractional-order parameter β has a stabilization effect.
(iv)
When s 2 = 0 and s 1 = 3 , 5 , 7 , this refers to the injection of external effector cells to destroy the tumor cells. From Figure 4, it follows that as the speed of tumor elimination becomes larger, the value of s 1 becomes larger. However, as the number increases, the change becomes less and less obvious. Which means that there is also a critical value when the value s 1 is selected.
(v)
When s 2 = 0 and s 1 = 3 , this indicates tumor therapy by external effector cells. When s 2 = 1 and s 1 = 3 , this indicates tumor therapy by external effector cells and external input of IL-2 into the model (Equation (1)). This shows that the injection of external IL-2 enhanced the speed of tumor extinction in Figure 5.
(vi)
In order to investigate the impact of external input of IL-2 into the model (Equation (1)) with different positive values, Figure 6 indicates that with a larger value of s 2 and the same value of s 1 , the speed of tumor extinction will be larger. Moreover, there is also a critical value when the value s 2 is selected.
(vii)
From Theorem 5, the tumor is incurable. However, immunotherapy can achieve an effective and stable tumor control. There is an equilibrium A ^ = ( 1.11 , 0.51 , 1.64 ) with a = 0.9 , m 1 = 1 , k 1 = 0.5 , b = 1 , c = 3 , s 1 = 0.3 , d = 5 , k 2 = 1 , m 2 = 1 , s 2 = 1.5 , and β = 0.9 . Figure 7 shows that the unique positive equilibrium A ^ is locally asymptotically stable.
(viii)
From Theorem 7, the tumor cells will be incurable. Figure 8 indicates that the unique positive equilibrium A ^ is always present and stable when the positive initial conditions are different; this shows that the tumor cells v cannot be permanently eradicated.
(ix)
In order to illustrate the expediency of using the fractal approach in comparison with the integer approach, Figure 9 shows that the integer approach is not practical, and the tumor cells are still growing with β = 2 , but the fractional approach is more realistic at this point and the number of tumor cells drops dramatically.

5. Conclusions

In this paper, we have investigated a fractional-order tumor–immune interaction model with immunotherapy and B-D functional response function. First, we obtain the existence, uniqueness, and nonnegativity of the solutions of the model (Equation (1)) (see Theorems 1 and 2). Second, there are some equilibrium points, depending on what treatment is not used and whether the tumor is curable. Moreover, the local asymptotic stability of these equilibrium points is proved by using the Hurwitz conditions and the linearization method (see Theorems 3–6), which means that the tumor cells can be destroyed or controlled by using immunotherapy. Then, sufficient conditions for the global asymptotic stability of the unique positive equilibrium are proved by using the method of the Lyapunov function (see Theorem 7). Finally, we visualized the theorem presented in this paper by numerical simulation. The numerical simulation results show that the equilibrium points can reach local or global asymptotic stability, which indicates that under certain conditions, the tumor can be completely cured or effectively controlled, respectively. For the fractional-order equation used in this paper, the larger the value of β , the faster and better the tumor will be cured or controlled. At present, the numerical simulation is still in the verification stage of theoretical results, and the simulation analysis has not involved real data. It is hoped that the real data can be better combined with theoretical research in the future, so as to achieve the effect of tumor suppression more efficiently.

Author Contributions

Conceptualization, X.F.; Investigation, X.F. and M.L.; Writing—original draft, M.L.; Writing—review and editing, X.F., Y.J. and D.L.; Visualization, X.F., M.L., Y.J. and D.L. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the National Natural Science Foundation of China grant number 12001425; the Social Science Key Research Plan in Shaanxi Province of China, China grant number 2022HZ1809; the National Innovation and Entrepreneurship Training Program for College Students, China grant number 202110702010; the Key Postgraduate Education Reform Project of Xi’an University of Technology, China grant number XAGDYJ220106.

Data Availability Statement

Not applicable.

Acknowledgments

The authors would like to thank the editor and reviewers greatly for their precious comments and suggestions.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Kawaguchi, K.; Maeshima, Y.; Toi, M. Tumor immune microenvironment and modelic response in breast cancer. Med. Oncol. 2022, 19, 1–10. [Google Scholar]
  2. Wu, Z.; Zhang, C.; Najafi, M. Targeting of the tumor immune microenvironment by metformin. J. Cell Commun. Signal. 2021, 16, 333–348. [Google Scholar] [CrossRef] [PubMed]
  3. Bashkirtseva, I.; Chukhareva, A.; Ryashko, L. Modeling and analysis of nonlinear tumor-immune interaction under chemotherapy and radiotherapy. Math. Methods Appl. Sci. 2021, 45, 7983–7991. [Google Scholar] [CrossRef]
  4. Lee, H.; Na, K.J.; Choi, H. Didderences in Tumor Immune Microenvironment in Matastatic Sites of Breast Cancer. Front. Oncol. 2021, 11, 722. [Google Scholar]
  5. Ko, W.; Ahn, I. Stationary patterns and stability in a tumor-immune interaction model with immunotherapy. J. Math. Anal. Appl. 2011, 383, 307–329. [Google Scholar] [CrossRef] [Green Version]
  6. Han, L.; Messan, M.R.; Yogurtcu, O.N.; Nukala, U.; Yang, H. Analysis of tumor-immune functional responses in a mathematical model of neoantigen cancer vaccines. Math. Biosci. 2023, 356, 108966. [Google Scholar] [CrossRef] [PubMed]
  7. Zhou, Q.; Yan, Y.; Li, Y.; Fu, H.; Lu, D.; Li, Z.; Wang, Y.; Wang, J.; Zhu, H.; Ren, J.; et al. Tumor-derived extracellular vesicles in melanoma immune response and immunotherapy. Biomed. Pharmacother. 2022, 156, 113790. [Google Scholar] [CrossRef] [PubMed]
  8. Kirschner, D.; Panetta, J.C. Modeling immunotherapy of the tumor—Immune interaction. J. Math. Biol. 1998, 37, 235–252. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  9. Yang, W.; Feng, X.; Liang, S.; Wang, X. Asymptotic Behavior Analysis of a Fractional-Order Tumor-Immune Interaction Model with Immunotherapy. Complexity 2020, 2020, 7062957. [Google Scholar] [CrossRef]
  10. Petráš, I. Fractional-Order Nonlinear Models: Modeling, Analysis and Simulation; Springer Science & Business Media: Berlin, Germany, 2011. [Google Scholar]
  11. Kilbas, A.A.A.; Srivastava, H.M.; Trujillo, J.J. Theory and Applications of Fractional Differential Equations; Elsevier Science Limited: Amsterdam, The Netherlands, 2006; Volume 204. [Google Scholar]
  12. Beddington, J.R. Mutual Interference Between Parasites or Predators and its Effect on Searching Efficiency. J. Anim. Ecol. 1975, 44, 331. [Google Scholar] [CrossRef] [Green Version]
  13. DeAngelis, D.; Goldstein, R.A.; Neill, R.V.O. A model for trophic interaction. Ecology 1975, 56, 881–892. [Google Scholar] [CrossRef]
  14. Li, H.-L.; Zhang, L.; Hu, C.; Jiang, Y.-L.; Teng, Z. Dynamical analysis of a fractional-order predator-prey model incorporating a prey refuge. J. Appl. Math. Comput. 2016, 54, 435–449. [Google Scholar] [CrossRef]
  15. Li, Y.; Chen, Y.; Podlubny, I. Stability of fractional-order nonlinear dynamic models: Lyapunov direct method and generalized Mittag-Leffler stability. Comput. Math. Appl. 2010, 59, 1810–1821. [Google Scholar] [CrossRef] [Green Version]
  16. Choi, S.K.; Kang, B.; Koo, N. Stability for caputo fractional differential model. Abstr. Appl. Anal. 2014, 2014, 631419. [Google Scholar] [CrossRef] [Green Version]
  17. Pao, C.V. Nonlinear Parabolic and Elliptic Equations; Plenum Press: New York, NY, USA, 1992. [Google Scholar]
  18. Clark, R.N. The Routh-Hurwitz stability criterion, revisited. IEEE Control Syst. 1992, 12, 119–120. [Google Scholar] [CrossRef] [Green Version]
  19. Vargas-De-León, C. Voletrra-type Lyapunov functions for fractional-order epidemic systems. Commun. Nonlinear Sci. Numer. Simul. 2015, 24, 75–85. [Google Scholar] [CrossRef]
  20. Aguila-Camacho, N.; Duarte-Mermoud, M.A.; Gallegos, J.A. Lyapunov functions for fractional order systems. Commun. Nonlinear Sci. Numer. Simul. 2014, 19, 2951–2957. [Google Scholar] [CrossRef]
  21. Diethelm, K.; Ford, N.J.; Freed, A.D. A Predictor-Corrector Approach for the Numerical Solution of Fractional Differential Equations. Nonlinear Dyn. 2002, 29, 3–22. [Google Scholar] [CrossRef]
  22. Garrappa, R. Trapezoidal methods for fractional differential equations: Theoretical and computational aspects. Math. Comput. Simul. 2013, 110, 96–112. [Google Scholar] [CrossRef] [Green Version]
Figure 1. Time series (a) and phase diagram (b) of equilibrium point A 2 = ( 3 , 0 , 0 ) with a = 0.9 , m 1 = 1 , k 1 = 0.5 , b = 1 , c = 3 , s 1 = 3 , d = 2.5 , k 2 = 1 , m 2 = 1 , s 2 = 0 , and β = 0.9 .
Figure 1. Time series (a) and phase diagram (b) of equilibrium point A 2 = ( 3 , 0 , 0 ) with a = 0.9 , m 1 = 1 , k 1 = 0.5 , b = 1 , c = 3 , s 1 = 3 , d = 2.5 , k 2 = 1 , m 2 = 1 , s 2 = 0 , and β = 0.9 .
Fractalfract 07 00200 g001
Figure 2. Time series (a) and phase diagram (b) of equilibrium point A 4 = ( 4 , 0 , 1 ) with a = 0.9 , m 1 = 1 , k 1 = 0.5 , b = 1 , c = 3 , s 1 = 3 , d = 2.5 , k 2 = 1 , m 2 = 1 , s 2 = 1 , and β = 0.9 .
Figure 2. Time series (a) and phase diagram (b) of equilibrium point A 4 = ( 4 , 0 , 1 ) with a = 0.9 , m 1 = 1 , k 1 = 0.5 , b = 1 , c = 3 , s 1 = 3 , d = 2.5 , k 2 = 1 , m 2 = 1 , s 2 = 1 , and β = 0.9 .
Fractalfract 07 00200 g002
Figure 3. Time series of model (Equation (1)) with a = 0.9 , m 1 = 1 , k 1 = 0.5 , b = 1 , c = 3 , d = 2.5 k 2 = 1 , m 2 = 1 , s 1 = 3 and different values of β . (a) s 2 = 0 and (b) s 2 = 1 .
Figure 3. Time series of model (Equation (1)) with a = 0.9 , m 1 = 1 , k 1 = 0.5 , b = 1 , c = 3 , d = 2.5 k 2 = 1 , m 2 = 1 , s 1 = 3 and different values of β . (a) s 2 = 0 and (b) s 2 = 1 .
Fractalfract 07 00200 g003
Figure 4. Time series of tumor cells with a = 0.9 , m 1 = 1 , k 1 = 0.5 , b = 1 , c = 3 , d = 2.5 , k 2 = 1 , m 2 = 1 , s 2 = 0 and different treatment by different densities of s 1 .
Figure 4. Time series of tumor cells with a = 0.9 , m 1 = 1 , k 1 = 0.5 , b = 1 , c = 3 , d = 2.5 , k 2 = 1 , m 2 = 1 , s 2 = 0 and different treatment by different densities of s 1 .
Fractalfract 07 00200 g004
Figure 5. Time series of tumor cells with a = 0.9 , m 1 = 1 , k 1 = 0.5 , b = 1 , c = 3 , d = 2.5 , k 2 = 1 , m 2 = 1 , s 1 = 3 and different treatment by the different densities of s 2 .
Figure 5. Time series of tumor cells with a = 0.9 , m 1 = 1 , k 1 = 0.5 , b = 1 , c = 3 , d = 2.5 , k 2 = 1 , m 2 = 1 , s 1 = 3 and different treatment by the different densities of s 2 .
Fractalfract 07 00200 g005
Figure 6. Time series of tumor cells with a = 0.9 , m 1 = 1 , k 1 = 0.5 , b = 1 , c = 3 , d = 2.5 , k 2 = 1 , m 2 = 1 , s 1 = 3 and different treatment by higher different densities of s 2 .
Figure 6. Time series of tumor cells with a = 0.9 , m 1 = 1 , k 1 = 0.5 , b = 1 , c = 3 , d = 2.5 , k 2 = 1 , m 2 = 1 , s 1 = 3 and different treatment by higher different densities of s 2 .
Fractalfract 07 00200 g006
Figure 7. Time series (a) and phase diagram (b) of equilibrium point A ^ = ( 1.11 , 0.51 , 1.64 ) with a = 0.9 , m 1 = 1 , k 1 = 0.5 , b = 1 , c = 3 , s 1 = 0.3 , d = 5 , k 2 = 1 , m 2 = 1 , s 2 = 1.5 , and β = 0.9 .
Figure 7. Time series (a) and phase diagram (b) of equilibrium point A ^ = ( 1.11 , 0.51 , 1.64 ) with a = 0.9 , m 1 = 1 , k 1 = 0.5 , b = 1 , c = 3 , s 1 = 0.3 , d = 5 , k 2 = 1 , m 2 = 1 , s 2 = 1.5 , and β = 0.9 .
Fractalfract 07 00200 g007
Figure 8. Phase diagram of equilibrium point A ^ = ( 1.11 , 0.51 , 1.64 ) with a = 0.9 , m 1 = 1 , k 1 = 0.5 , b = 1 , c =   3 , s 1 = 0.3 , d = 5 , k 2 = 1 , m 2 = 1 , s 2 = 1.5 , β = 0.9 , and different initial conditions ( u 0 , v 0 , w 0 ) . The initial conditions of the blue, red and green lines are ( u 0 , v 0 , w 0 ) = ( 4 , 4 , 4 ) , ( u 0 , v 0 , w 0 ) = ( 3 , 3 , 3 ) and ( u 0 , v 0 , w 0 ) = ( 2 , 2 , 2 ) respectively.
Figure 8. Phase diagram of equilibrium point A ^ = ( 1.11 , 0.51 , 1.64 ) with a = 0.9 , m 1 = 1 , k 1 = 0.5 , b = 1 , c =   3 , s 1 = 0.3 , d = 5 , k 2 = 1 , m 2 = 1 , s 2 = 1.5 , β = 0.9 , and different initial conditions ( u 0 , v 0 , w 0 ) . The initial conditions of the blue, red and green lines are ( u 0 , v 0 , w 0 ) = ( 4 , 4 , 4 ) , ( u 0 , v 0 , w 0 ) = ( 3 , 3 , 3 ) and ( u 0 , v 0 , w 0 ) = ( 2 , 2 , 2 ) respectively.
Fractalfract 07 00200 g008
Figure 9. Time series of tumor cells with a = 0.9 , m 1 = 1 , k 1 = 0.5 , b = 1 , c = 3 , d = 2.5 , k 2 = 1 , m 2 = 1 , s 1 = 3 , s 2 = 1 and different treatment by higher different densities of β .
Figure 9. Time series of tumor cells with a = 0.9 , m 1 = 1 , k 1 = 0.5 , b = 1 , c = 3 , d = 2.5 , k 2 = 1 , m 2 = 1 , s 1 = 3 , s 2 = 1 and different treatment by higher different densities of β .
Fractalfract 07 00200 g009
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Feng, X.; Liu, M.; Jiang, Y.; Li, D. Dynamics and Stability of a Fractional-Order Tumor–Immune Interaction Model with B-D Functional Response and Immunotherapy. Fractal Fract. 2023, 7, 200. https://doi.org/10.3390/fractalfract7020200

AMA Style

Feng X, Liu M, Jiang Y, Li D. Dynamics and Stability of a Fractional-Order Tumor–Immune Interaction Model with B-D Functional Response and Immunotherapy. Fractal and Fractional. 2023; 7(2):200. https://doi.org/10.3390/fractalfract7020200

Chicago/Turabian Style

Feng, Xiaozhou, Mengyan Liu, Yaolin Jiang, and Dongping Li. 2023. "Dynamics and Stability of a Fractional-Order Tumor–Immune Interaction Model with B-D Functional Response and Immunotherapy" Fractal and Fractional 7, no. 2: 200. https://doi.org/10.3390/fractalfract7020200

APA Style

Feng, X., Liu, M., Jiang, Y., & Li, D. (2023). Dynamics and Stability of a Fractional-Order Tumor–Immune Interaction Model with B-D Functional Response and Immunotherapy. Fractal and Fractional, 7(2), 200. https://doi.org/10.3390/fractalfract7020200

Article Metrics

Back to TopTop