Time-fractional diffusion model on dynamical effect of dendritic cells on HIV pathogenesis

Abstract

Chronic immune activation and viral latency are among the main factors behind the persistence of Human Immunodeficiency Virus (HIV) and other viral diseases. Immune activation is chiefly lead by dendritic cells which detect the presence of Virus (HIV) and notify other immune cells particularly the helper T cells via the process called antigen presentation. However, antigen presentation process increases the chances of T cells to be infected by the viral particles trapped by dendritic cells. In this paper, a time-fractional diffusion model is proposed to investigate the impact of dendritic cells on the spread of HIV among susceptible T cells with respect to time and anomalous diffusion posed by cells crowding. The equilibrium points of the model are obtained and analysed. The disease free steady state proved to be stable both in spatially homogeneous case and in the presence diffusion and chemotaxis. The endemic steady state is stable in the absence of diffusion and chemotaxis, however the presence of diffusion induces instability if a thresshold value is exceeded. A priori estimates were also obtained in the appropriate Sobolev spaces. Furthermore, the numerical experiments were conducted to examine the dynamic behaviour of cells densities with respect to time and the sub-diffusion parameter.

Keywords

Dendritic cells time fractional model HIV antigen presentation

1. Introduction

Dendritic cells (DCs) play multiple roles at different stages of HIV infection characterized by microbe sensing and communication with T cell (antigen presentation) via soluble and cell-associated molecules. Only antigen-presenting cells such as macrophages, B lymphocytes, and DCs can activate a naive helper T cell when the matching antigen is presented [17, 21], because every helper T cell is dedicated to a specific antigen. However, macrophages and B cells can only activate memory T cells whereas DCs can activate both memory and naive T cells. These make DCs the most potent of all the antigen-presenting cells [17].

Unfortunately, their very close interactions with CD4+ T cells create an avenue for the infection of new target cells. The HIV virions are captured by DCs using C-C chemokine receptor type 5 (CCR5) and C-X-C chemokine receptor type 4 (CXCR4) and retain infectious viruses for long periods [21]. There are strong evidences that DCs are not productively infected by HIV; rather, DCs serve as a reservoir of virions at their surface with high probability of infecting the target CD4+ T cells [47]. DC maturation is associated with a diminished ability to support HIV replication as compared to immature DC [18].

The mystery behind persistent progress of HIV infection lies in chronic immune activation lead by DCs and viral reservoirs. However, research on HIV infection particularly amongst computational scientists is majorly directed to examining key variables responsible for the progress of the disease. For instance, mathematical modellers were initially interested to study the dynamics of the susceptible target cells, infected target cells and the free viral particles [22, 52]. Even though the model is able to fit data, major factors such as immune response and other transmission mechanisms were not included. In the last four decades, there had been growing body of literature dedicated to extended models that capture virus-to-cell transmission [57, 50], direct cell-to-cell transmission [34, 56], antiretroviral therapies [39], latency of infected CD4+ T cells [8], drug resistance [44, 48] and time delays [38, 51].

Mathematical modelling of immune response to HIV is still an underdeveloped research area with few articles found in literature which capture CD8+ responses [27, 35, 53], the interactions between virus, CD4+ T cells and CD8+ T cells [54]. All these models focus only on time evolution while the spatial distribution of interacting cells receives less or no attention. In addition, the impact of DCs on progression and pathology of HIV cannot be overemphasised and it remains a gap in mathematical modelling of infectious disease. Truly, war against most viral infections such as HIV cannot be won without a clear understanding of immune cells activation and mobilization.

Fractional calculus is now an indispensable modelling tool in science and engineering systems that deviate from the classical laws [32]. Examples of such models are found in hydrology [12], plasma turbulence [10], finance [7], epidemiology [6, 20, 14, 45] and ecology [9, 19]. Indeed, the modelling of HIV epidemic using fractional derivatives has received reasonable attention over the past two decades. Mainly, these time fractional models focus on basic target cells dynamics [16], viral latency [41], therapy [3] and time delay [59] with respect to time evolution.

The diffusion of a particle in a simple medium is modelled using classical Fick’s law where mean-square displacement of a diffusing particle is linearly scaled by $\left\langle r^{2}\left(t\right)\right\rangle\sim t$ . Whereas, in complex media such as lymph nodes, blood streams and other tissues, anomalous diffusion is expected where the mean-square displacement obeys the power law $\langle r^{2}(t)\rangle\sim t^{\alpha}$ . A diffusion is called subdiffusion if $0<\alpha<1$ , superdiffusion when $1<\alpha<2$ and normal when $\alpha=1$ . In the review [23], crowding in biological media is identified as one of the major causes of subdiffusion in biological systems along side with random traps, strong interaction of particles in single-file, etc. [4, 29, 28]. In this paper, we propose a time-fractional diffusion model to study the effects of antigen presentation by DCs on the progress of HIV, capturing subdiffusion in time and chemotactic movement of DCs.

2. Preliminaries

In this section we will recall some basic definitions and results that are useful in the sequel.

.

(Fractional Integrals) [43]: Let $I=[a,b]\mathrm{}\subset\mathbb{R}$ be a finite interval. The left and right Riemann-Liouville fractional integrals of order $\alpha\in{\mathbb{R}}^{+}$ of a function $v(t)$ denoted by ${}_{a}{J^{\alpha}_{t}}v(t)$ and ${}_{t}{J^{\alpha}_{b}}v(t)$ , are defined by

${}_{a}{J^{\alpha}_{t}}v(t)=\frac{1}{\Gamma(\alpha)}\int^{t}_{a}{{(t-\tau)}^{% \alpha-1}v(\tau)d\tau}∼{}∼{}\alpha>0,t>a,$ (1) ${}_{t}{J^{\alpha}_{b}}v(t)=\frac{1}{\Gamma(\alpha)}\int^{b}_{t}{{(\tau-t)}^{% \alpha-1}v(\tau)d\tau}∼{}∼{}\alpha>0,b>t,$ (2)

.

(Fractional Derivatives) [43]: Caputo and Riemann-Liouville fractional left and right hands derivatives of order $\alpha\in{\mathbb{R}}^{+},$ $n-1<\alpha<n,$ $n\in\mathbb{N}$ of a function $v\left(x,t\right)$ are given as

${}^{C}_{a}{{\partial}^{\alpha}_{t}}v(x,t)=\frac{1}{\Gamma(n-\alpha)}\int^{t}_{% a}{\frac{{\partial}^{n}v(x,t)}{\partial{\tau}^{n}}{(t-\tau)}^{n-\alpha-1}d\tau% },∼{}∼{}\alpha>0,t>a,$ (3) ${}^{C}_{t}{{\partial}^{\alpha}_{b}}v(x,t)=\frac{(-1)^{n}}{\Gamma(n-\alpha)}% \int^{b}_{t}{\frac{{\partial}^{n}v(x,t)}{\partial{\tau}^{n}}{(\tau-t)}^{n-% \alpha-1}d\tau},∼{}∼{}\alpha>0,b>t,$ (4)

and

${}^{R}_{a}{{\partial}^{\alpha}_{t}}v(x,t)=\frac{1}{\Gamma(n-\alpha)}\frac{{% \partial}^{n}}{\partial t^{n}}\int^{t}_{a}{v(x,\tau){(t-\tau)}^{n-\alpha-1}d% \tau},∼{}∼{}\alpha>0,t>a,$ (5) ${}^{R}_{t}{{\partial}^{\alpha}_{b}}v(x,t)=\frac{(-1)^{n}}{\Gamma(n-\alpha)}% \frac{{\partial}^{n}}{\partial t^{n}}\int^{b}_{t}{v(x,\tau){(\tau-t)}^{n-% \alpha-1}d\tau},∼{}∼{}\alpha>0,b>t,$ (6)

Note: In the paper, we denote ${}^{C}_{0}{{\partial}^{\alpha}_{t}}={}^{C}{{\partial}^{\alpha}_{t}}$ and ${}^{R}_{0}{{\partial}^{\alpha}_{t}}={}^{R}{{\partial}^{\alpha}_{t}}.$

.

For an integer $m>0$ and real $p$ with $1\leqslant p<\infty$ , the Sobolev space is defined by [1]

$W^{m,p}(\Omega)=\{v\in L^{p}(\Omega)|D^{\alpha}v\in L^{p}(\Omega)\forall|% \alpha|\leqslant m\}\$ (7)

equipped with the following norms

$\displaystyle{\left\|v\right\|}_{W^{m,p}\left(\Omega\right)}\coloneqq\left\{% \begin{array}[]{l}{\left(\sum_{\left|\alpha\right|\leqslant m}{\int_{\Omega}{{% \left|D^{\alpha}v\right|}^{p}dx}}\right)}^{{1}/{p}}1\leqslant p<\infty\\ \sum_{\left|\alpha\right|\leqslant m}{\textit{ess}{\mathop{\sup}_{\Omega}\left% |D^{\alpha}v\right|}}∼{}∼{}p=\infty\end{array}\right.$ (8) $\displaystyle{\left|v\right|}_{W^{m,p}\left(\Omega\right)}\coloneqq{\left(\sum% _{\left|\alpha\right|=m}{\int_{\Omega}{{\left|D^{\alpha}v\right|}^{p}dx}}% \right)}^{{1}/{p}}1\leqslant p<\infty.$ (9)

Now, for $p=2$ , a Hilbert space is defined $W^{m,2}(\Omega)=H^{m}(\Omega)$ with the inner product

${(u,v)}_{m,\Omega}=\sum_{|\alpha|\leqslant m}{{(D^{\alpha}u,D^{\alpha}v)}_{0,% \Omega}},\quad H^{1}_{0}(\Omega)=\{v\in H^{1}|v=0\text{ on }\partial\Omega\}.$ (10)

.

The one and two parameters Mittag-Lefler functions are given as [30]

$\displaystyle E_{\alpha}(v)=\sum^{\infty}_{k=0}{\frac{v^{k}}{\Gamma(\alpha k+1% )}},∼{}∼{}\alpha>0,$ (11) $\displaystyle E_{\alpha,\beta}(v)=\sum^{\infty}_{k=0}{\frac{v^{k}}{\Gamma(% \alpha k+\beta)}},∼{}∼{}\alpha>0,\beta>0,$ (12)

.

Useful Inequalities

Young’s inequality: Let $p, q$ be positive real numbers satisfying $\frac{1}{p}+\frac{1}{q}=1$ then,

$\displaystyle ab\leqslant\frac{a^{p}}{p}+\frac{b^{q}}{q},\text{ for }a,b% \geqslant 0,$ $\displaystyle ab\leqslant\frac{a^{2}}{2}+\frac{b^{2}}{2},\text{ for }p=q=2.$ (13)

Poincaré inequality: Let $1\leqslant p<\infty$ and $\Omega$ a subset with at least one bound. Then there exists a constant $C(p,\Omega)$ , such that, for every function $v\in H^{1}_{0}(\Omega)$

${\|v\|}_{L^{p}(\Omega)}\leqslant C{\|\nabla v\|}_{L^{p}(\Omega)}.$ (14)

Generalized Grönwall’s inequality [58]: Suppose $\alpha>0$ , $x(t)$ is a nonnegative function locally integrable on $0\leqslant t<T$ (some $T\leqslant+\infty$ ) and $y(t)$ is a nonnegative, nondecreasing continuous function defined on $0\leqslant t<T$ , $y(t)\leqslant M$ (constant), and suppose $v(t)$ is nonnegative and locally integrable on $0\leqslant t<T$ with

$v(t)\leqslant x(t)+y(t)\int^{t}_{0}{{(t-s)}^{\alpha-1}v(s)ds},∼{}∼{}0\leqslant t% <T,$ (15)

then

$v(t)\leqslant x(t)+y(t)\int^{t}_{0}{\left[\sum^{\infty}_{k=1}{\frac{{(y(t)% \Gamma(\alpha))}^{k}}{\Gamma(k\alpha)}{(t-s)}^{k\alpha-1}x(s)}\right]ds},∼{}∼{% }0\leqslant t<T.$ (16)

3. Mathematical model

To account for the effect of anomalous diffusion (subdiffusion) of immune cells due to cell crowding, we will capture the nonlinear time scaling phenomenon (with $0<\alpha<1$ ) modelled by the left side Caputo fractional derivative defined by Eq. (3) above. Certainly, the mean-square displacement $\langle r^{2}(t)\rangle\sim t^{\alpha}$ of the interacting cells may differ. But to minimize the complexity of the model associated with immune response and subdiffusion, we assume ${\alpha}$ common to all cells. The Caputo derivative is preferred over Riemann-Liouville derivative because Riemann-Liouville derivative cannot be coupled with classical initial conditions of the form we considered in this model [13]. Before now, several investigators have successfully used Caputo derivative for epidemiological models [11, 16, 14, 45].

In this model, the unknowns include: $c_{1}(x,t)$ : Average population naive immature DCs; $c_{2}(x,t)$ : Average population immature DCs with bound viral particle for presentation; $c_{3}(x,t)$ : Average population mature DCs; $u_{1}(x,t)$ : Average population of naive Helper T cells; $u_{2}(x,t)$ : Average population of actively infected Helper T cells; $v(x,t)$ : Average population density of HIV.

3.1 Dendritic cells

DCs are called matured only when they complete the antigen presentation process. Naive immature DCs move randomly in the body to keep watch against invasion of pathogens. Once a viral particle is located, the DC binds and takes it to helper T cells. Here, we modelled in Eq. (17) ${\tilde{D}}_{1}\Delta c_{1}$ and ${\rho}_{1}$ as the diffusion and source term for naive immature DCs respectively, ${\tilde{b}}_{3}c_{3}$ represent a number of mature DCs that revert to naïve state after antigen presentation process. ${\tilde{b}}_{1}c_{1}v$ and ${\widetilde{\gamma}}_{1}c_{1}$ are naive immature DCs activated due to contact with virons and death of naive immature DCs respectively. In Eq. (18), ${\tilde{D}}_{2}\Delta c_{2}$ is diffusion of partially mature DCs; $\widetilde{\chi}\nabla\cdot(c_{2}\nabla u_{1})$ is chemotaxis of partially mature DCs towards naive helper T cells; ${\tilde{b}}_{1}c_{1}v$ is a portion of population of naive immature DCs that are partially matured; $({\tilde{b}}_{2}+{\widetilde{\gamma}}_{2})c_{2}$ is the sum of partially mature DCs that underwent antigen presentation process and dead of partially mature DCs. Equation (19) captures the diffusion of mature DCs as ${\tilde{D}}_{3}\Delta c_{3}$ ; ${\tilde{b}}_{2}c_{2}$ is proportion of partially mature DCs that moves into fully matured state; $({\tilde{b}}_{3}+{\widetilde{\gamma}}_{3})c_{3}$ is the sum of fully matured DCs that revert to naive state and death of mature DCs. Hence the equations are given with initial-boundary conditions,

$\displaystyle{}^{C}{{\partial}^{\alpha}_{t}}c_{1}={\tilde{D}}_{1}\Delta c_{1}+% {\rho}_{1}+{\tilde{b}}_{3}c_{3}-{\tilde{b}}_{1}c_{1}v-{\widetilde{\gamma}}_{1}% c_{1},$ (17) $\displaystyle{}^{C}{{\partial}^{\alpha}_{t}}c_{2}={\tilde{D}}_{2}\Delta c_{2}-% \widetilde{\chi}\nabla\cdot(c_{2}\nabla u_{1})+{\tilde{b}}_{1}c_{1}v-({\tilde{% b}}_{2}+{\widetilde{\gamma}}_{2})c_{2},$ (18) $\displaystyle{}^{C}{{\partial}^{\alpha}_{t}}c_{3}={\tilde{D}}_{3}\Delta c_{3}+% {\tilde{b}}_{2}c_{2}-({\tilde{b}}_{3}+{\widetilde{\gamma}}_{3})c_{3},$ (19) $\displaystyle c_{i}(x,0)=c^{0}_{i}(x),\frac{\partial c_{i}}{\partial n}=0,∼{}∼% {}i=1,2,3.$

3.2 T-cells

In this model, the T-cells are divided into susceptible ( $u_{1}$ ) and infected ( $u_{2}$ ). In Eq. (20), ${\tilde{D}}_{4}\Delta u_{1}$ is the diffusion of susceptible T-cells and ${\rho}_{2}$ represents the source terms. Respectively, ${\tilde{k}}_{1}u_{1}v$ and ${\widetilde{\mu}}_{1}u_{1}$ are numbers of infected and dead susceptible T-cells. The term ${\tilde{k}}_{3}c_{2}u_{1}$ captures the number of susceptible T cells infected via contact with DCs during antigen presentation. In the same manner, in Eq. (21), ${\tilde{D}}_{5}\Delta u_{2}$ is the diffusion term, ${\tilde{k}}_{1}u_{1}v$ and ${\tilde{k}}_{3}c_{2}u_{1}$ are infection terms, $({\tilde{k}}_{2}+{\widetilde{\mu}}_{2})u_{2}$ accounts for lysing and death terms respectively.

$\displaystyle{}^{C}{{\partial}^{\alpha}_{t}}u_{1}={\tilde{D}}_{4}\Delta u_{1}+% {\rho}_{2}-{\tilde{k}}_{1}u_{1}v-{\tilde{k}}_{3}c_{2}u_{1}-{\widetilde{\mu}}_{% 1}u_{1},$ (20) $\displaystyle{}^{C}{{\partial}^{\alpha}_{t}}u_{2}={\tilde{D}}_{5}\Delta u_{2}+% {\tilde{k}}_{1}u_{1}v+{\tilde{k}}_{3}c_{2}u_{1}-({\tilde{k}}_{2}+{\widetilde{% \mu}}_{2})u_{2},$ (21) $\displaystyle u_{i}(x,0)=u^{0}_{i}(x),\frac{\partial u_{i}}{\partial n}=0,∼{}∼% {}i=1,2.$

3.3 Virus

Here, ${\tilde{D}}_{6}\Delta v$ represents the diffusion of the viral particles; $N{\tilde{k}}_{2}u_{2}$ is the replication of virus via lysing of actively infected T cells. The virus clearance rate is ${\widetilde{\mu}}_{3}$ .

${}^{C}{{\partial}^{\alpha}_{t}}v={\tilde{D}}_{6}\Delta v+N{\tilde{k}}_{2}u_{2}% -{\widetilde{\mu}}_{3}v,$

(22) $\displaystyle v(x,0)=v^{0}(x),\frac{\partial v}{\partial n}=0.$

Note that, the normal diffusion system is recovered for $\alpha=1$ . However, system Eqs (17)–(3.3) posed dimension inconsistency in the time variable. That is, the left hand side of the system has dimension ${(\textit{time})}^{-\alpha}$ while the right hand side has dimension ${(\textit{time})}^{-1}$ . This is mathematically corrected by setting

$\displaystyle b_{1}={\tilde{b}}^{\alpha}_{1},b_{2}={\tilde{b}}^{\alpha}_{2},b_% {3}={\tilde{b}}^{\alpha}_{3},k_{1}={\tilde{k}}^{\alpha}_{1},k_{2}={\tilde{k}}^% {\alpha}_{2},{\gamma}_{1}={\widetilde{\gamma}}^{\alpha}_{1},{\gamma}_{2}={% \widetilde{\gamma}}^{\alpha}_{2},{\gamma}_{3}={\widetilde{\gamma}}^{\alpha}_{3% },{\mu}_{1}={\widetilde{\mu}}^{\alpha}_{1},{\mu}_{2}={\widetilde{\mu}}^{\alpha% }_{2},$ $\displaystyle{\mu}_{3}={\widetilde{\mu}}^{\alpha}_{3},D_{1}={\tilde{D}}^{% \alpha}_{1},D_{2}={\tilde{D}}^{\alpha}_{2},D_{3}={\tilde{D}}^{\alpha}_{3},D_{4% }={\tilde{D}}^{\alpha}_{4},D_{5}={\tilde{D}}^{\alpha}_{5},D_{6}={\tilde{D}}^{% \alpha}_{6},\chi={\widetilde{\chi}}^{\alpha}.$

This correction technique was introduced by [14]. Hence, the correct form of the model is given as

${}^{C}{{\partial}^{\alpha}_{t}}c_{1}=D_{1}\Delta c_{1}+{\rho}_{1}+b_{3}c_{3}-{% b_{1}c}_{1}v-{\gamma}_{1}c_{1},$ (23) ${}^{C}{{\partial}^{\alpha}_{t}}c_{2}=D_{2}\Delta c_{2}-\chi\nabla\cdot(c_{2}% \nabla u_{1})+{b_{1}c}_{1}v-(b_{2}+{\gamma}_{2})c_{2},$ (24) ${}^{C}{{\partial}^{\alpha}_{t}}c_{3}=D_{3}\Delta c_{3}+b_{2}c_{2}-(b_{3}+{% \gamma}_{3})c_{3},$ (25) ${}^{C}{{\partial}^{\alpha}_{t}}u_{1}=D_{4}\Delta u_{1}+{\rho}_{2}-k_{3}c_{2}u_% {1}-k_{1}u_{1}v-{\mu}_{1}u_{1},$ (26) ${}^{C}{{\partial}^{\alpha}_{t}}u_{2}=D_{5}\Delta u_{2}+k_{3}c_{2}u_{1}+k_{1}u_% {1}v-(k_{2}+{\mu}_{2})u_{2},$ (27) ${}^{C}{{\partial}^{\alpha}_{t}}v=D_{6}\Delta v+Nk_{2}u_{2}-{\mu}_{3}v,$ (28) $\displaystyle c_{i}(x,0)=c^{0}_{i}(x),\frac{\partial c_{i}}{\partial n}=0,∼{}∼% {}i=1,2,3,$ $\displaystyle u_{i}(x,0)=u^{0}_{i}(x),\frac{\partial u_{i}}{\partial n}=0,∼{}∼% {}i=1,2,$ $\displaystyle v(x,0)=v^{0}(x),\frac{\partial v}{\partial n}=0.$

where $\frac{\partial c_{i}}{\partial n},\frac{\partial u_{i}}{\partial n},\frac{% \partial v}{\partial n}$ represents the zero-flux boundary conditions. The model parameters are described in Table 1 above.

Table 1
Parameter values

Symbol Description Values & references

${\tilde{b}}_{1}$ Rate of activation of DCs 0.3/day [25]

${\tilde{b}}_{2}$ Rate of maturation of DCs 0.4/day [25]

${\tilde{b}}_{3}$ Rate which matured DCs revert to naive state 0.04/day

${\tilde{k}}_{1}$ Free virus-to-cell infection rate 0.05 virons mm ${}^{3}$ /day [55]

${\tilde{k}}_{2}$ Rate of lysing actively infected Helper T cells 0.25/day [55]

${\tilde{k}}_{3}$ Infection rate via contact with DCs 0.1/day

${\rho}_{1}$ Source term for naive immature DCs $100$ cells[24]

${\rho}_{2}$ Source term for susceptible Helper T cells $100$ cells[24]

${\widetilde{\gamma}}_{1}$ Death rate of immature DCs 0.35/day [26]

${\widetilde{\gamma}}_{2}$ Death rate of partially mature DCs 0.25/day [26]

${\widetilde{\gamma}}_{3}$ Death rate of fully mature DCs 1.0/day [26]

${\widetilde{\mu}}_{1}$ Death rate of susceptible Helper T cells 0.06/day

${\widetilde{\mu}}_{2}$ Death rate of infected Helper T cells 0.24/day

${\widetilde{\mu}}_{3}$ Death rate of free HIV 3/day

${\tilde{D}}_{1}$ Diffusion rate of DCs 0.0028 mm ${}^{2}$ /day [31]

${\tilde{D}}_{2}$ Diffusion rate of partially mature DCs 0.0028 mm ${}^{2}$ /day [31]

${\tilde{D}}_{3}$ Diffusion rate of fully mature DCs 0.0028 mm ${}^{2}$ /day [33]

$\chi$ Chemotaxis rate of partially mature DCs 0.0028 mm ${}^{2}$ /day [33]

${\tilde{D}}_{4}$ Diffusion rate of susceptible helper T cells 0.0045 mm ${}^{2}$ /day [37]

${\tilde{D}}_{5}$ Diffusion rate of infected helper T cells 0.0045 mm ${}^{2}$ /day[37]

${\tilde{D}}_{6}$ Diffusion rate of HIV 0.0008 mm ${}^{2}$ /day

$N$ Number of viruses produced per infected helper T cells 10 virons/cell [24]

4. Analysis

Symbol	Description	Values & references
${\tilde{b}}_{1}$	Rate of activation of DCs	0.3/day [25]
${\tilde{b}}_{2}$	Rate of maturation of DCs	0.4/day [25]
${\tilde{b}}_{3}$	Rate which matured DCs revert to naive state	0.04/day
${\tilde{k}}_{1}$	Free virus-to-cell infection rate	0.05 virons mm ${}^{3}$ /day [55]
${\tilde{k}}_{2}$	Rate of lysing actively infected Helper T cells	0.25/day [55]
${\tilde{k}}_{3}$	Infection rate via contact with DCs	0.1/day
${\rho}_{1}$	Source term for naive immature DCs	$100$ cells[24]
${\rho}_{2}$	Source term for susceptible Helper T cells	$100$ cells[24]
${\widetilde{\gamma}}_{1}$	Death rate of immature DCs	0.35/day [26]
${\widetilde{\gamma}}_{2}$	Death rate of partially mature DCs	0.25/day [26]
${\widetilde{\gamma}}_{3}$	Death rate of fully mature DCs	1.0/day [26]
${\widetilde{\mu}}_{1}$	Death rate of susceptible Helper T cells	0.06/day
${\widetilde{\mu}}_{2}$	Death rate of infected Helper T cells	0.24/day
${\widetilde{\mu}}_{3}$	Death rate of free HIV	3/day
${\tilde{D}}_{1}$	Diffusion rate of DCs	0.0028 mm ${}^{2}$ /day [31]
${\tilde{D}}_{2}$	Diffusion rate of partially mature DCs	0.0028 mm ${}^{2}$ /day [31]
${\tilde{D}}_{3}$	Diffusion rate of fully mature DCs	0.0028 mm ${}^{2}$ /day [33]
$\chi$	Chemotaxis rate of partially mature DCs	0.0028 mm ${}^{2}$ /day [33]
${\tilde{D}}_{4}$	Diffusion rate of susceptible helper T cells	0.0045 mm ${}^{2}$ /day [37]
${\tilde{D}}_{5}$	Diffusion rate of infected helper T cells	0.0045 mm ${}^{2}$ /day[37]
${\tilde{D}}_{6}$	Diffusion rate of HIV	0.0008 mm ${}^{2}$ /day
$N$	Number of viruses produced per infected helper T cells	10 virons/cell [24]

4.1 Equilibrium states

By setting the kinetic parts of the equation to zero and solving simultaneously, the disease free and endemic equilibrium states of the system Eqs (23)–(28) are obtained respectively:

$\displaystyle E_{1}[c^{*}_{1},c^{*}_{2},c^{*}_{3},u^{*}_{1},u^{*}_{2},v^{*}]=% \left(\frac{{\rho}_{1}}{{\gamma}_{1}},0,0,\frac{{\rho}_{2}}{{\mu}_{1}},0,0% \right),$ (29) $\displaystyle E_{2}[c^{*}_{1},c^{*}_{2},c^{*}_{3},u^{*}_{1},u^{*}_{2},v^{*}]=% \left(\frac{\rho_{1}}{(b_{1}v^{*}-\gamma_{1})}+\frac{b_{1}b_{2}b_{3}\rho_{1}}{% B},\frac{b_{1}\rho_{1}(b_{3}+\gamma_{3})}{B},\frac{b_{1}b_{2}\rho_{1}}{B},% \frac{1}{M},\right.$ $\displaystyle∼{}\left.\mspace{190.0mu }\frac{1}{\mu_{3}(k_{2}+\mu_{2})},\frac{% Nk_{2}u_{2}}{\mu_{3}}\right),$ (30)

where $B=(b_{3}+\gamma_{3})(b_{2}+\gamma_{2})(b_{1}v^{*}-\gamma_{1})-b_{1}b_{2}b_{3},$ and $M=\frac{b_{1}\rho_{1}\mu_{3}k_{3}(b_{3}+\gamma_{3})}{B}+\frac{Nk_{2}k_{1}}{\mu% _{3}(k_{2}+\mu_{2})}$ . We proceed with the stability analysis of equilibrium points in the following theorems.

.

The disease free equilibrium point $E_{1}[c^{*}_{1},c^{*}_{2},c^{*}_{3},u^{*}_{1},u^{*}_{2},v^{*}]$ of the system Eqs (23)–(28) is uniformly and asymptotically stable in case of the kinetic system and remains stable in the presence of diffusion.

Proof..

The absence of the virus reduces system into

${}^{C}{{\partial}^{\alpha}_{t}}c_{1}=D_{1}\Delta c_{1}+{\rho}_{1}-{\gamma}_{1}% c_{1},\quad{}^{C}{{\partial}^{\alpha}_{t}}c_{2}=0,\quad{}^{C}{{\partial}^{% \alpha}_{t}}c_{3}=0,$ (31) ${}^{C}{{\partial}^{\alpha}_{t}}u_{1}=D_{4}\Delta u_{1}+{\rho}_{2}-{\mu}_{1}u_{% 1},\quad{}^{C}{{\partial}^{\alpha}_{t}}u_{1}=0,\quad{}^{C}{{\partial}^{\alpha}% _{t}}v=0,$ (32)

which leads to the homogenous disease free equilibrium point $E_{1}=\left(\frac{{\rho}_{1}}{{\gamma}_{1}},0,0,\frac{{\rho}_{2}}{{\mu}_{1}},0% ,0\right).$ Simple calculations give leads to two negative eigenvalues, i.e.

$\left|\begin{array}[]{c@{\quad}c}-{\gamma}_{1}-\lambda\quad&0\\ 0\quad&-{\mu}_{1}-\lambda\end{array}\right|=0,∼{}∼{}\Longrightarrow{\lambda}_{% 1}=-{\gamma}_{1},∼{}∼{}{\lambda}_{2}=-{\mu}_{1}.$ (33)

This implies that the kinetic system of Eqs (31)–(32) is stable. To examine diffusion induced instability of the system, Eqs (31)–(32) is rewritten as

$\displaystyle 0=D_{1}\Delta c_{1}+\rho_{1}-{\gamma}_{1}c_{1},$ (34) $\displaystyle 0=D_{4}\Delta u_{1}+{\rho}_{2}-{\mu}_{1}u_{1}.$ (35)

Let ${\varphi}_{n}={\cos({\pi nx}/{l})}{\sin({\pi nx}/{l})}$ be the eigenvector of the Laplacian together with the eigenvalues $-k^{2}_{n}=-{({\pi nx}/{l})}^{2}$ on a bounded square. Linearization about $E_{1}$ will reduce to $2\times 2$ Jacobian matrices

$J\left(k\right)=\begin{pmatrix}-k^{2}_{n}D_{1}-{\gamma}_{1}&0\\ 0&-k^{2}_{n}D_{4}-{\mu}_{1}\end{pmatrix}.$ (36)

This matrix has the following eigenvalues

${\lambda}_{1}=-k^{2}_{n}D_{1}-{\gamma}_{1},∼{}∼{}{\lambda}_{2}=-k^{2}_{n}D_{4}% -{\mu}_{1},\quad\text{ for all }k>0.\$ (37)

Since all the eigenvalues are negative, the system remains stable in the presence of diffusion. ∎

.

(Gershgorin circle theorem) [5]: Let $A$ be a $n\times n$ matrix, with entries $a_{ij}$ . For $i=1,\dots,n$ , let $R_{i}=\sum_{i\neq j}{|a_{ij}|}$ be the sum of the absolute values of the non-diagonal entries in the i-th row and that $D(a_{ij},R_{i})$ be the closed disc centered at $a_{ij}$ with radius $R_{i}$ . Then every eigenvalue of $A$ lies within at least one of the Gershgorin discs $D(a_{ij},R_{i}).$

.

The endemic equilibrium point $E_{2}[c^{*}_{1},c^{*}_{2},c^{*}_{3},u^{*}_{1},u^{*}_{2},v^{*}]$ of the system is uniformly and asymptotically stable for $c^{*}_{1},c^{*}_{2},c^{*}_{3},u^{*}_{1},u^{*}_{2},v^{*}\geqslant 0.$ However, if $\frac{b_{1}(c^{*}_{1}+v^{*})+k^{2}_{n}\chi c^{*}_{2}}{k^{2}_{n}D_{2}}>1$ the diffusion-chemotaxis system becomes unstable with the possibility of Turing patterns.

Proof..

Assume that the diffusion and the chemotactic movement of cells are negligible. Then we obtain the Jacobian matrix of the system Eqs (23)–(28) as

$J=\begin{pmatrix}-b_{1}v^{*}-{\gamma}_{1}&0&b_{3}&0&0&-b_{1}c^{*}_{1}\\ b_{1}v^{*}&-b_{2}-{\gamma}_{2}&0&0&0&b_{1}c^{*}_{1}\\ 0&b_{2}&-\left(b_{3}+{\gamma}_{3}\right)&0&0&0\\ 0&-k_{3}u^{*}_{1}&0&-k_{3}c^{*}_{2}-k_{1}v^{*}-{\mu}_{1}&0&-k_{1}u^{*}_{1}\\ 0&k_{3}u^{*}_{1}&0&k_{3}c^{*}_{2}+k_{1}v^{*}&-\left(k_{2}+{\mu}_{2}\right)&k_{% 1}u^{*}_{1}\\ 0&0&0&0&Nk_{2}&-{\mu}_{3}\end{pmatrix}.$ (38)

We need to prove that the real parts of the eigenvalues of are strictly negative. Using parameter values in Table 1, for $E_{2}[c^{*}_{1},c^{*}_{2},c^{*}_{3},u^{*}_{1},u^{*}_{2},v^{*}]$ we compute the eigenvalues numerically by setting

$|J(E_{2})-\lambda I|=0$ (39)

and solving characteristic polynomial, we obtained

$\lambda_{1}=-4.1;{\lambda}_{2}=-3.4;{\lambda}_{3}=-1.6;{\lambda}_{4}=-1.1;{% \lambda}_{5}=-0.62;{\lambda}_{6}=-0.31.$ (40)

Since all the eigenvalues are negative, we deduce that the endemic equilibrium point of the system is uniformly and asymptotically stable. Next, the linearization of the full diffusion-chemotaxis time independent system of Eqs (23)–(28) leads to matrix $A$ below which is the sum of $A_{0}$ and $J$ pertaining to diffusion-chemotaxis and kinetic terms respectively.

$A(k^{2})=\begin{pmatrix}-k^{2}_{n}D_{1}-a_{11}&0&b_{3}&0&0&-b_{1}c^{*}_{1}\\ b_{1}v^{*}&-k^{2}_{n}D_{2}-a_{22}&0&k^{2}_{n}\chi c^{*}_{2}&0&b_{1}c^{*}_{1}\\ 0&b_{2}&-k^{2}_{n}D_{3}-a_{33}&0&0&0\\ 0&-k_{3}u^{*}_{1}&0&-k^{2}_{n}D_{4}-a_{44}&0&-k_{1}u^{*}_{1}\\ 0&k_{3}u^{*}_{1}&0&k_{3}c^{*}_{2}+k_{1}v^{*}&-k^{2}_{n}D_{5}-a_{55}&k_{1}u^{*}% _{1}\\ 0&0&0&0&Nk_{2}&-k^{2}_{n}D_{6}-a_{66}\end{pmatrix}$ (41)

where $J$ is given by Eq. (38) and $A_{0}$ is given as

$A_{0}(k^{2})=\begin{pmatrix}-k^{2}_{n}D_{1}&0&0&0&0&0\\ 0&-k^{2}_{n}D_{2}&0&\chi c^{*}_{2}&0&0\\ 0&b_{2}&-k^{2}_{n}D_{3}&0&0&0\\ 0&0&0&-k^{2}_{n}D_{4}&0&0\\ 0&0&0&0&-k^{2}_{n}D_{5}&0\\ 0&0&0&0&0&-k^{2}_{n}D_{6}\end{pmatrix}.$ (42)

The stability of the system is examined using Gershgorin circle theorem. We evaluate the Gershgorin disc for rows $i=1,\dots,6$ as follows

$\displaystyle i=1:|{\lambda}_{1}+k^{2}_{n}D_{1}+a_{11}|<b_{3}+b_{1}c^{*}_{1},$ (43) $\displaystyle i=2:|{\lambda}_{2}+k^{2}_{n}D_{2}+a_{22}|<b_{1}v^{*}+k^{2}_{n}% \chi c^{*}_{2}+b_{1}c^{*}_{1},$ (44) $\displaystyle i=3:|{\lambda}_{3}+k^{2}_{n}D_{3}+a_{33}|<b_{2},$ (45) $\displaystyle i=4:|{\lambda}_{4}+k^{2}_{n}D_{4}+a_{44}|<k_{1}u^{*}_{1}+k_{3}u^% {*}_{1},$ (46) $\displaystyle i=5:|{\lambda}_{5}+k^{2}_{n}D_{5}+a_{55}|<k_{3}c^{*}_{2}+k_{3}u^% {*}_{1}+k_{1}v^{*}+k_{1}u^{*}_{1},$ (47) $\displaystyle i=6:|{\lambda}_{6}+k^{2}_{n}D_{6}+a_{66}|<Nk_{2},$ (48)

where $a_{11}=b_{1}v^{*}+{\gamma}_{1},a_{22}=b_{2}+{\gamma}_{2},a_{33}=b_{3}+{\gamma}% _{3},a_{44}=k_{3}c^{*}_{2}+k_{1}v^{*}+{\mu}_{1},a_{55}=k_{2}+{\mu}_{2},a_{66}=% {\mu}_{3}.$

All the eigenvalues are centered on the negative part of the real axis. However, by numerical inspection using the parameter values in Table 1, we discover that the Gershgorin disc corresponding to row $i=2$ contains 0 in the interior. This may lead to diffusion induced instability with depending on activation, diffusion and chemotaxis rates of partially mature DCs and $k^{2}_{n}$ with the following condition

$\frac{b_{1}(c^{*}_{1}+v^{*})+k^{2}_{n}\chi c^{*}_{2}}{k^{2}_{n}D_{2}}>1.$ (49)

With this, we conclude that the stable endemic equilibrium point $E_{2}[c^{*}_{1},c^{*}_{2},c^{*}_{3},u^{*}_{1},u^{*}_{2},v^{*}]$ turns out to be unstable in presence of diffusion and chemotaxis provided Eq. (49) holds, hence the proof. ∎

4.2 A priori estimates

.

Let $\Omega\mathrm{}$ be a bounded domain in ${\mathbb{R}}^{2}$ with smooth boundary. Given that $c^{0}_{1}>0,c^{0}_{2}\geqslant 0,c^{0}_{3}\geqslant 0,u^{0}_{1}>0,u^{0}_{2}% \geqslant,v^{0}>0$ in $\overline{\Omega}$ , $(c^{0}_{1},c^{0}_{2},c^{0}_{3},u^{0}_{1},u^{0}_{2},v^{0})\in{(H^{1}(\Omega))}^% {6}={\bm{H}}^{1}_{0}(\Omega)$ , and its dual ${\bm{H}}^{-1}(\Omega)$ , then system Eqs (23)–(28) admit a unique solution for $0\leqslant t\leqslant T$

$\bm{U}\in\bm{W}(0,T)\coloneqq\{\bm{U}\in{\bm{L}}^{2}(0,T;{\bm{H}}^{1}_{0}(% \Omega));{}^{C}{{\partial}^{\alpha}_{t}}\bm{U}\in{\bm{L}}^{2}(0,T;{\bm{H}}^{-1% }(\Omega))\},\$ (50)

where $\bm{U}=\left(c_{1},c_{2},c_{3},u_{1},u_{2},v\right)$ .

In addition, for constants $\epsilon,{\lambda}_{0},\xi,K>0,$

$\displaystyle\mathop{\sup}_{t\in[0,T]}\psi(t)\leqslant\epsilon E_{\alpha}({% \lambda}_{0}T^{\alpha})+K,$ (51) $\displaystyle{\mathop{\sup}_{t\in[0,T]}\varphi(t)}\leqslant\frac{\varphi(0)}{1% -\frac{\xi T^{\alpha}}{\Gamma(\alpha+1)}},$ (52)

where $E_{\alpha}\left(t\right)$ is one parameter Mittag-Lefler function and

$\displaystyle\psi(t)\coloneqq({\|c_{1}(t)\|}^{2}_{L^{2}}+{\|c_{2}(t)\|}^{2}_{L% ^{2}}+{\|c_{2}(t)\|}^{2}_{L^{2}}+{\|u_{1}(t)\|}^{2}_{L^{2}}+{\|u_{2}(t)\|}^{2}% _{L^{2}}+{\|v(t)\|}^{2}_{L^{2}}),$ (53) $\displaystyle\varphi(t)\coloneqq({\|\nabla c_{1}(t)\|}^{2}_{L^{2}}+{\|\nabla c% _{2}(t)\|}^{2}_{L^{2}}+{\|\nabla c_{3}(t)\|}^{2}_{L^{2}}+{\|\nabla u_{1}(t)\|}% ^{2}_{L^{2}}+{\|\nabla u_{2}(t)\|}^{2}_{L^{2}}+{\|\nabla v(t)\|}^{2}_{L^{2}}).$ (54)

Proof..

The main goal of the proof is to obtain a priori estimates with which boundedness and uniqueness will be checked.

Estimate I: Multiplying Eq. (17) by $c_{1}$ and integrating over the domain $\Omega$ , we have

$\frac{1}{2}{}^{C}{{\partial}^{\alpha}_{t}}\int_{\Omega}{c^{2}_{1}}=D_{1}\int_{% \Omega}{c_{1}\Delta c_{1}}+{\rho}_{1}\int_{\Omega}{c_{1}}+b_{3}\int_{\Omega}{c% _{1}c_{3}}-b_{1}\int_{\Omega}{vc^{2}_{1}}-{\gamma}_{1}\int_{\Omega}{c^{2}_{1}},$ (55)

integrating the diffusion term by parts and applying the boundary conditions

$\frac{1}{2}{}^{C}{{\partial}^{\alpha}_{t}}\int_{\Omega}{c^{2}_{1}}+D_{1}\int_{% \Omega}{{(\nabla c_{1})}^{2}}\leqslant{\rho}_{1}\int_{\Omega}{c^{2}_{1}}+b_{3}% \int_{\Omega}{c_{1}c_{3}}-b_{1}\int_{\Omega}{vc^{2}_{1}}-{\gamma}_{1}\int_{% \Omega}{c^{2}_{1}}.$ (56)

Applying the young’s inequality,

$\frac{1}{2}{}^{C}{{\partial}^{\alpha}_{t}}\int_{\Omega}{c^{2}_{1}}+D_{1}\int_{% \Omega}{{(\nabla c_{1})}^{2}}\leqslant({\rho}_{1}-{\gamma}_{1})\int_{\Omega}{c% ^{2}_{1}}+\frac{b_{3}}{2}\int_{\Omega}{(c^{2}_{1}+c^{2}_{3})}-b_{1}\int_{% \Omega}{vc^{2}_{1}},$ (57)

since $b_{1}\geqslant 0$ and ${\rho}_{1}>{\gamma}_{1},$

${}^{C}{{\partial}^{\alpha}_{t}}{\|c_{1}(t)\|}^{2}_{L^{2}}+2D_{1}{\|\nabla c_{1% }(t)\|}^{2}_{L^{2}}\leqslant 2({\rho}_{1}+b_{3}-{\gamma}_{1}){\|c_{1}(t)\|}^{2% }_{L^{2}}+2b_{3}{\|c_{3}(t)\|}^{2}_{L^{2}}.$ (58)

Secondly, multiplying Eq. (18) by $c_{2}$ and integrating over the domain $\Omega$ with boundary conditions in mind, we obtain

$\int_{\Omega}{c_{2}{}^{C}{{\partial}^{\alpha}_{t}}c_{2}}=D_{2}\int_{\Omega}{c_% {2}\Delta c_{2}}-{\chi}_{1}\int_{\Omega}{c_{2}\nabla(c_{2}\nabla u_{1})}+b_{1}% \int_{\Omega}{vc_{1}c_{2}}-(b_{2}+{\gamma}_{2})\int_{\Omega}{c_{2}c_{2}}.$ (59)

The diffusion and chemotaxis terms are integrated by parts

$\frac{1}{2}{}^{C}{{\partial}^{\alpha}_{t}}\int_{\Omega}{c^{2}_{2}}+D_{2}\int_{% \Omega}{{(\nabla c_{2})}^{2}}\leqslant{\chi}_{1}\int_{\Omega}{c_{2}\nabla u_{1% }.\nabla c_{2}}+b_{1}\int_{\Omega}{vc_{1}c_{2}}-(b_{2}+{\gamma}_{2})\int_{% \Omega}{c^{2}_{2}},$ (60)

but by Young’s inequality

$\int_{\Omega}{c_{2}\nabla u_{1}.\nabla c_{2}}\leqslant\int_{\Omega}{\left(% \frac{{(c_{2}\nabla u_{1})}^{2}}{2}+\frac{{(\nabla c_{2})}^{2}}{2}\right)}% \leqslant\int_{\Omega}{\left(\frac{C_{1}c^{2}_{2}}{2}+\frac{{(\nabla c_{2})}^{% 2}}{2}\right)},$ (61)

which follow from Sobolev embedding for $\nabla u_{1}\in L^{2}(\Omega)$ , the $|\nabla u_{1}|\leqslant C_{1}$ .

Similarly, using the Young’s inequality, it leads to

$\int_{\Omega}{vc_{1}c_{2}}\leqslant\int_{\Omega}{\left(\frac{v^{2}c^{2}_{1}}{2% }+\frac{c^{2}_{2}}{2}\right)}.$

By Cauchy-Schwarz inequality this leads to

${}^{C}{{\partial}^{\alpha}_{t}}{\|c_{2}(t)\|}^{2}_{L^{2}}+2\left(D_{2}-\frac{{% \chi}_{1}}{2}\right){\|\nabla c_{2}(t)\|}^{2}_{L^{2}}\leqslant\left(b_{1}+% \frac{C_{1}}{2}-b_{2}-{\gamma}_{2}\right){\|c_{2}(t)\|}^{2}_{L^{2}}+$ $\displaystyle\quad∼{}b_{1}{\|v(t)\|}^{2}_{L^{2}}{\|c_{1}(t)\|}^{2}_{L^{2}}.$ (62)

Next, we follow the same argument for the remaining equations to have

${}^{C}{{\partial}^{\alpha}_{t}}{\|c_{3}(t)\|}^{2}_{L^{2}}\leqslant 2D_{3}{\|% \nabla c_{3}(t)\|}^{2}_{L^{2}}+b_{2}{\|c_{2}(t)\|}^{2}_{L^{2}}+(b_{2}-b_{3}-{% \gamma}_{3}){\|c_{3}(t)\|}^{2}_{L^{2}},$ (63) ${}^{C}{{\partial}^{\alpha}_{t}}{\|u_{1}(t)\|}^{2}_{L^{2}}\leqslant 2D_{4}{\|% \nabla u_{1}(t)\|}^{2}_{L^{2}}+2({\rho}_{2}-{\mu}_{1}){\|u_{1}(t)\|}^{2}_{L^{2% }}+k_{3}{\|c_{2}(t)\|}^{2}_{L^{2}},$ (64) ${}^{C}{{\partial}^{\alpha}_{t}}{\|u_{2}(t)\|}^{2}_{L^{2}}+2D_{5}{\|\nabla u_{2% }(t)\|}^{2}_{L^{2}}\leqslant 2\left(\frac{k_{1}}{2}-k_{2}-{\textrm{\`{i}}}_{2}% \right){\|u_{2}(t)\|}^{2}_{L^{2}}+$ $\displaystyle\quad∼{}k_{1}{\|v(t)\|}^{2}_{L^{2}}{\|u_{1}(t)\|}^{2}_{L^{2}}+k_{% 3}{\|c_{2}(t)\|}^{2}_{L^{2}},$ (65) ${}^{C}{{\partial}^{\alpha}_{t}}{\|v(t)\|}^{2}_{L^{2}}+2D_{6}{\|\nabla v(t)\|}^% {2}_{L^{2}}\leqslant 2(Nk_{2}-{\mu}_{3}){\|v(t)\|}^{2}_{L^{2}}+Nk_{2}{\|u_{2}(% t)\|}^{2}_{L^{2}}.$ (66)

Now, adding Eqs (58), (4.2) and (66), ignoring the negative terms, we arrive at

${}^{C}{{\partial}^{\alpha}_{t}}({\left\|c_{1}\left(t\right)\right\|}^{2}_{L^{2% }}+{\left\|c_{2}\left(t\right)\right\|}^{2}_{L^{2}}+{\left\|c_{2}\left(t\right% )\right\|}^{2}_{L^{2}}+{\left\|u_{1}\left(t\right)\right\|}^{2}_{L^{2}}+{\left% \|u_{2}\left(t\right)\right\|}^{2}_{L^{2}}+{\left\|v\left(t\right)\right\|}^{2% }_{L^{2}})$ $\displaystyle\leqslant\lambda_{0}({\left\|c_{1}\left(t\right)\right\|}^{2}_{L^% {2}}+{\left\|c_{2}\left(t\right)\right\|}^{2}_{L^{2}}+{\left\|c_{3}\left(t% \right)\right\|}^{2}_{L^{2}}+{\left\|u_{1}\left(t\right)\right\|}^{2}_{L^{2}}+% +{\left\|u_{2}\left(t\right)\right\|}^{2}_{L^{2}}+{\left\|v\left(t\right)% \right\|}^{2}_{L^{2}})+{\lambda}_{1},$ (67) $\displaystyle{\lambda}_{0}\coloneqq{\max\left\{2\left({\rho}_{1}+b_{3}\right),% \left(b_{1}+b_{2}+\frac{C_{1}}{2}\right),\left(b_{2}+2b_{3}\right),2{\rho}_{2}% ,Nk_{2}+k_{1},2Nk_{2}\right\}},$ $\displaystyle{\lambda}_{1}\coloneqq b_{1}{\left\|v^{0}\right\|}^{2}_{L^{2}}{% \left\|c^{0}_{1}\right\|}^{2}_{L^{2}}+k_{1}{\left\|v^{0}\right\|}^{2}_{L^{2}}{% \left\|u^{0}_{1}\right\|}^{2}_{L^{2}},$

where we used the fact that

${\|v(t)\|}^{2}_{L^{2}}\leqslant C{\|v^{0}\|}^{2}_{L^{2}},{\|c_{1}(t)\|}^{2}_{L% ^{2}}\leqslant C{\|c^{0}_{1}\|}^{2}_{L^{2}},{\|u_{1}(t)\|}^{2}_{L^{2}}% \leqslant C{\|u^{0}_{1}\|}^{2}_{L^{2}}.$

Then, using the substitution Eq. (53), we take the fractional integral on both sides of Eq. (67) and apply generalized Grönwall’s lemma to realize

${\mathop{\sup}_{t\in[0,T]}\psi(t)}\leqslant\epsilon E_{\alpha}({\lambda}_{0}T^% {\alpha})+K,$ (68)

where the constants $\epsilon=\epsilon(c^{0}_{1},c^{0}_{2},c^{0}_{3},u^{0}_{1},u^{0}_{2},v^{0},% \alpha,T)$ while $K=K({\lambda}_{1},\alpha,T)$ and

Estimate II: Multiplying Eq. (17) by $-\Delta c_{1}$ and integrating over the domain $\Omega$ implies

$-\frac{1}{2}{}^{C}{{\partial}^{\alpha}_{t}}\int_{\Omega}{c_{1}\Delta c_{1}}=-D% _{1}\int_{\Omega}{\Delta c_{1}\cdot\Delta c_{1}}-{\rho}_{1}\int_{\Omega}{% \Delta c_{1}}-b_{3}\int_{\Omega}{c_{3}\Delta c_{1}}+b_{1}\int_{\Omega}{vc_{1}% \Delta c_{1}}+{\gamma}_{1}\int_{\Omega}{c_{1}\Delta c_{1}}.$ (69)

Integrating all the terms by parts except the diffusion term and applying the boundary conditions we have

$\displaystyle\frac{1}{2}{}^{C}{{\partial}^{\alpha}_{t}}\int_{\Omega}{\nabla c_% {1}\cdot\nabla c_{1}}=-D_{1}\int_{\Omega}{\Delta c_{1}.\Delta c_{1}}-{\rho}_{1% }\int_{\Omega}{\Delta c_{1}}+b_{3}\int_{\Omega}{\nabla c_{3}\cdot\nabla c_{1}}% -b_{1}$ $\displaystyle\int_{\Omega}{\nabla(vc_{1})\cdot\nabla c_{1}}-{\gamma}_{1}\int_{% \Omega}{\nabla c_{1}.\nabla c_{1}},$ (70) ${}^{C}{{\partial}^{\alpha}_{t}}{\|\nabla c_{1}(t)\|}^{2}_{L^{2}}+2(D_{1}+{\rho% }_{1}){\|\Delta c_{1}(t)\|}^{2}_{L^{2}}\leqslant(b_{3}-{\gamma}_{1}){\|\nabla c% _{1}(t)\|}^{2}_{L^{2}}+b_{3}{\|\nabla c_{3}(t)\|}^{2}_{L^{2}},$ (71)

where we dropped the negative term and the Young’s inequality is applied to have

$\int_{\Omega}{\nabla c_{3}\cdot\nabla c_{1}}\leqslant\int_{\Omega}{\left(\frac% {{(\nabla c_{3})}^{2}}{2}\mathrm{+}\frac{(\nabla c_{1})^{2}}{2}\right)}.$

Secondly, we multiplying Eq. (18) by $-\Delta c_{2}$ and integrating over the domain $\Omega$

$\displaystyle-\frac{1}{2}{}^{C}{{\partial}^{\alpha}_{t}}\int_{\Omega}{\Delta c% _{2}c_{2}}=-D_{2}\int_{\Omega}{\Delta c_{2}\cdot\Delta c_{2}}+{\chi}_{1}\int_{% \Omega}{\Delta c_{2}\cdot\nabla(c_{2}\nabla u_{1})}-$ $\displaystyle\mspace{157.0mu }b_{1}\int_{\Omega}{vc_{1}\Delta c_{2}}+(b_{2}+{% \gamma}_{2})\int_{\Omega}{c_{2}\Delta c_{2}},$ (72) $\displaystyle\frac{1}{2}{}^{C}{{\partial}^{\alpha}_{t}}\int_{\Omega}{\nabla c_% {2}\cdot\nabla c_{2}}+D_{2}\int_{\Omega}{\Delta c_{2}\cdot\Delta c_{2}}={\chi}% _{1}\int_{\Omega}{\Delta c_{2}\cdot\nabla(c_{2}\nabla u_{1})}+b_{1}\int_{% \Omega}{\nabla(vc_{1})\cdot\nabla c_{2}}-$ $\displaystyle\mspace{314.0mu }(b_{2}+{\gamma}_{2})\int_{\Omega}{\nabla c_{2}% \cdot\nabla c_{2}},$ (73)

but it is easy to see that

$\int_{\Omega}{\nabla(vc_{1})\cdot\nabla c_{2}}=\int_{\Omega}{(c_{1}\nabla v+v% \nabla c_{1})\cdot\nabla c_{2}}=\int_{\Omega}{(c_{1}\nabla v\cdot\nabla c_{2}+% v\nabla c_{1}\cdot\nabla c_{2})}.$ (74)

Due to the fact that, for $\nabla u_{1},\nabla v,\nabla c_{1}\in L^{2}(\Omega)$ , then $|\nabla u_{1}|\leqslant C_{1},|\nabla v|\leqslant C_{2}$ and $|\nabla c_{1}|\leqslant C_{3}$ , we conclude that

${}^{C}{{\partial}^{\alpha}_{t}}{\|\nabla c_{2}(t)\|}^{2}_{L^{2}}+\left(2D_{1}-% \frac{{\chi}_{1}}{2}\right){\|\Delta c_{2}(t)\|}^{2}_{L^{2}}\leqslant(b_{1}-b_% {2}-{\gamma}_{2}){\|\nabla c_{2}(t)\|}^{2}_{L^{2}}+$ $\displaystyle\quad∼{}b_{1}C_{3}{\|v(t)\|}^{2}_{L^{2}}+b_{1}C_{2}{\|c_{1}(t)\|}% ^{2}_{L^{2}}.$ (75)

Note that, the Young’s inequality was again invoked to disassemble the products.

Now, we proceed with the same series of argument to obtain estimate for the remaining equation however, details are omitted. This yields the following equations:

${}^{C}{{\partial}^{\alpha}_{t}}{\|\nabla c_{3}(t)\|}^{2}_{L^{2}}+2D_{3}{\|% \Delta c_{3}(t)\|}^{2}_{L^{2}}\leqslant b_{2}{\|\nabla c_{2}(t)\|}^{2}_{L^{2}}% +(b_{2}-b_{3}-{\gamma}_{3}){\|\nabla c_{3}(t)\|}^{2}_{L^{2}},$ (76) ${}^{C}{{\partial}^{\alpha}_{t}}{\|\nabla u_{1}(t)\|}^{2}_{L^{2}}+2(D_{4}+{\rho% }_{2}){\|\Delta u_{1}(t)\|}^{2}_{L^{2}}\leqslant{\mu}_{1}{\|\nabla u_{1}(t)\|}% ^{2}_{L^{2}}+k_{3}{\|\nabla c_{2}(t)\|}^{2}_{L^{2}},$ (77) ${}^{C}{{\partial}^{\alpha}_{t}}{\|\nabla u_{2}(t)\|}^{2}_{L^{2}}+2D_{5}{\|% \Delta u_{2}(t)\|}^{2}_{L^{2}}\leqslant k_{1}{\|\nabla u_{2}(t)\|}^{2}_{L^{2}}% +k_{1}C_{3}{\|v(t)\|}^{2}_{L^{2}}+$ $\displaystyle\quad∼{}(k_{1}C_{2}-k_{2}-{\mu}_{2}){\|u_{1}(t)\|}^{2}_{L^{2}},$ (78) ${}^{C}{{\partial}^{\alpha}_{t}}{\|\nabla v(t)\|}^{2}_{L^{2}}+2D_{6}{\|\Delta v% (t)\|}^{2}_{L^{2}}\leqslant Nk_{2}{\|\nabla u_{2}(t)\|}^{2}_{L^{2}}+(Nk_{2}-{% \mu}_{3}){\|\nabla v(t)\|}^{2}_{L^{2}}.$ (79)

Combining Eqs (71), (4.2)–(79) and dropping the negative terms

${}^{C}{{\partial}^{\alpha}_{t}}({\|\nabla c_{1}(t)\|}^{2}_{L^{2}}+{\|\nabla c_% {2}(t)\|}^{2}_{L^{2}}+{\|\nabla c_{3}(t)\|}^{2}_{L^{2}}+{\|\nabla u_{1}(t)\|}^% {2}_{L^{2}}+{\|\nabla u_{2}(t)\|}^{2}_{L^{2}}+{\|\nabla v(t)\|}^{2}_{L^{2}})$ $\displaystyle\leqslant{\xi}_{0}({\|\nabla c_{1}(t)\|}^{2}_{L^{2}}+{\|\nabla c_% {2}(t)\|}^{2}_{L^{2}}+{\|\nabla c_{3}(t)\|}^{2}_{L^{2}}+{\|\nabla u_{1}(t)\|}^% {2}_{L^{2}}+{\|\nabla u_{2}(t)\|}^{2}_{L^{2}}+{\|\nabla v(t)\|}^{2}_{L^{2}})$ $\displaystyle\quad+∼{}\xi_{1}({\|c_{1}(t)\|}^{2}_{L^{2}}+{\|u_{1}(t)\|}^{2}_{L% ^{2}}+{\|v(t)\|}^{2}_{L^{2}}).$ (80)

Thanks to Poincaré inequality, for a constant ${\xi}_{2}$

$({\|c_{1}(t)\|}^{2}_{L^{2}}+{\|u_{1}(t)\|}^{2}_{L^{2}}+{\|v(t)\|}^{2}_{L^{2}})% \leqslant{\xi}_{2}({\|\nabla c_{1}(t)\|}^{2}_{L^{2}}+{\|\nabla u_{1}(t)\|}^{2}% _{L^{2}}+{\|\nabla v(t)\|}^{2}_{L^{2}}),$ (81)

so that the right hand side of Eq. (81) can easily be absorbed in Eq. (80) which leads to

${}^{C}{{\partial}^{\alpha}_{t}}({\|\nabla c_{1}(t)\|}^{2}_{L^{2}}+{\|\nabla c_% {2}(t)\|}^{2}_{L^{2}}+{\|\nabla c_{3}(t)\|}^{2}_{L^{2}}+{\|\nabla u_{1}(t)\|}^% {2}_{L^{2}}+{\|\nabla u_{2}(t)\|}^{2}_{L^{2}}+{\|\nabla v(t)\|}^{2}_{L^{2}})$ $\displaystyle\leqslant\xi({\|\nabla c_{1}(t)\|}^{2}_{L^{2}}+{\|\nabla c_{2}(t)% \|}^{2}_{L^{2}}+{\|\nabla c_{3}(t)\|}^{2}_{L^{2}}+{\|\nabla u_{1}(t)\|}^{2}_{L% ^{2}}+{\|\nabla u_{2}(t)\|}^{2}_{L^{2}}{+}{\|\nabla v(t)\|}^{2}_{L^{2}}),$ (82)

where

${\xi}_{0}\coloneqq{\max(b_{3},b_{1}+b_{2},b_{2}+b_{3},{\mu}_{1},Nk_{2}+k_{1},% Nk_{2})},{\xi}_{1}\coloneqq\max(b_{1}C_{2},k_{1}C_{3}),\xi={\xi}_{0}+{\xi}_{1}.$

Using the substitution in Eq. (54), Eq. (82) is rewritten as

${}^{C}{{\partial}^{\alpha}_{t}}\varphi(t)\leqslant\xi\varphi(t).$ (83)

Taking the fractional integral on both sides of Eq. (83), we have

$\displaystyle\varphi(t)-\varphi(0)\leqslant\frac{\xi}{\Gamma(\alpha)}\int^{t}_% {0}{{(t-s)}^{\alpha-1}\varphi(s)ds},$ (84) $\displaystyle{\mathop{\sup}_{t\in[0,T]}\varphi(t)}\leqslant\varphi(0)+\frac{% \xi T^{\alpha}}{\Gamma(\alpha+1)}{\mathop{\sup}_{t\in[0,T]}\varphi(t)},$ (85) $\displaystyle\therefore{\mathop{\sup}_{t\in[0,T]}\varphi(t)}\leqslant\frac{% \varphi(0)}{1-\frac{\xi T^{\alpha}}{\Gamma(\alpha+1)}}.$ (86)

Uniqueness: Uniqueness of the solution of the system can be investigated by assuming that the system has two solutions $({\tilde{c}}_{1},{\tilde{c}}_{2},{\tilde{c}}_{3},{\tilde{u}}_{1},{\tilde{u}}_{% 2},\tilde{v})$ and $({\overline{c}}_{1},{\overline{c}}_{2},{\overline{c}}_{3},{\overline{u}}_{1},{% \overline{u}}_{2},\overline{v})$ with the same initial conditions $({\tilde{c}}^{0}_{1},{\tilde{c}}^{0}_{2},{\tilde{c}}^{0}_{3},{\tilde{u}}^{0}_{% 1},{\tilde{u}}^{0}_{2},{\tilde{v}}^{0})\in{\bm{H}}^{1}(\Omega)$ . Let $C_{1}={\tilde{c}}_{1}-{\overline{c}}_{1},C_{2}={\tilde{c}}_{2}-{\overline{c}}_% {2},C_{3}={\tilde{c}}_{3}-{\overline{c}}_{3},U_{1}={\tilde{u}}_{1}-{\overline{% u}}_{1},$ $U_{2}={\tilde{u}}_{2}-{\overline{u}}_{2},V=\tilde{v}-\overline{v}$ with $(C_{1},C_{2},C_{3},U_{1},U_{2},V)\in\bm{W}(0,T)$ satisfying Eqs (23)–(28). Now, using the estimate Eq. (67) we have

${}^{C}{{\partial}^{\alpha}_{t}}({\|C_{1}(t)\|}^{2}_{L^{2}}+{\|C_{2}(t)\|}^{2}_% {L^{2}}+{\|C_{2}(t)\|}^{2}_{L^{2}}+{\|U_{1}(t)\|}^{2}_{L^{2}}+{\|U_{2}(t)\|}^{% 2}_{L^{2}}+{\|V(t)\|}^{2}_{L^{2}})$ $\displaystyle\leqslant{\lambda}_{0}({\|C_{1}(t)\|}^{2}_{L^{2}}+{\|C_{2}(t)\|}^% {2}_{L^{2}}+{\|C_{3}(t)\|}^{2}_{L^{2}}+{\|U_{1}(t)\|}^{2}_{L^{2}}+{\|U_{2}(t)% \|}^{2}_{L^{2}}+{\|V(t)\|}^{2}_{L^{2}})+{\lambda}_{1},$ (87) $\displaystyle{\lambda}_{0}\coloneqq{\max\left\{2({\rho}_{1}+b_{3}),\left(b_{1}% +b_{2}+\frac{K_{1}}{2}\right),(b_{2}+2b_{3}),2{\rho}_{2},Nk_{2}+k_{1},2Nk_{2}% \right\}},$ $\displaystyle{\lambda}_{1}\coloneqq b_{1}{\|{\tilde{v}}^{0}\|}^{2}_{L^{2}}{\|{% \tilde{c}}^{0}_{1}\|}^{2}_{L^{2}}+k_{1}{\|{\tilde{v}}^{0}\|}^{2}_{L^{2}}{\|{% \tilde{u}}^{0}_{1}\|}^{2}_{L^{2}}.$

The desired result follows upon the application of generalized Grönwall’s lemma, hence the proof. ∎

5. Numerical scheme

In this section, we used the Grunwald-Letnikov formula in time [40], central in space scheme [42] to solve the system. Fractional derivatives hold feature of nonlocality that makes them difficult to discretize. In literature, Caputo derivative is often discretized through the definition of Grunwald-Letnikov fractional derivative [40, 46, 15] and via quadrature approximations of the fractional integral [60, 2]. This can also be achieved by generalized Taylor series expansion where fractional Euler scheme is obtained [49, 36]. Now, we take advantage of the relationship between Grunwald-Letnikov, Riemann-Liouville and Caputo fractional derivatives to obtain the desired scheme. Let $t_{n}=t_{0}+n\tau,$ where $t_{0}$ is the initial time and $\tau$ is the step size, then Grunwald-Letnikov fractional derivative is given by

${}^{GL}{{\partial}^{\alpha}_{t}}v(x,t)={\mathop{\lim}_{\tau\downarrow 0}\frac{% 1}{\tau}\sum^{\left[\frac{t-t_{0}}{\tau}\right]}_{n=0}{{(-1)}^{n}\binom{\alpha% }{n}}}v(x,t-n\tau)$ (88)

where $\binom{\alpha}{n}=\frac{\Gamma(\alpha+1)}{\Gamma(n+1)\Gamma(n-\alpha+1)}$ . For $0<\alpha<1$ , we recall the relationship between Riemann-Liouville and Caputo fractional derivatives

${}^{C}{{\partial}^{\alpha}_{t}}v(x,t)={}^{R}{{\partial}^{\alpha}_{t}}v(x,t)-% \frac{v(x,0)t^{-\alpha}}{\Gamma(1-\alpha)}.$ (89)

Using an equidistant grid in $[t_{0},T]$ , we have Grunwald-Letnikov approximation given by

$\frac{1}{{\tau}^{\alpha}}{\Delta}^{\alpha}_{\tau}v(x,t)=\frac{1}{{\tau}^{% \alpha}}\sum^{k+1}_{n=0}{{\beta}^{\alpha}_{n}v(x,t_{k+1-n})}$ (90)

where ${\beta}^{\alpha}_{n}$ are the fractional binomial coefficients $\binom{\alpha}{n}$ with the recursion formula

$\beta^{\alpha}_{n}=\left(1-\frac{1+\alpha}{n}\right){\beta}^{\alpha}_{n-1},∼{}% ∼{}{\beta}^{\alpha}_{0}=1.$ (91)

Now, the Riemann-Liouville fractional derivative is defined in terms Grunwald-Letnikov approximation formula as

${}^{R}{{\partial}^{\alpha}_{t}}v(x,t)={\mathop{\mathrm{lim}}_{\tau\to 0}\frac{% 1}{{\tau}^{\alpha}}{\Delta}^{\alpha}_{\tau}v(x,t)}.$ (92)

In view of Eqs (90)–(92), the Caputo fractional derivative is approximated as

${}^{C}{{\partial}^{\alpha}_{t}}v(x,t_{n+1})\approx\frac{1}{{\tau}^{\alpha}}% \sum^{k+1}_{n=0}{{\beta}^{\alpha}_{n}v(x,t_{k+1-n})}-{\tilde{v}}_{n+1},$ (93) $\displaystyle{\tilde{v}}_{n+1}=\frac{v(x,0)t^{-\alpha}_{n+1}}{\Gamma(1-\alpha)}.$ (94)

The correction term Eq. (94) vanishes with homogeneous initial conditions. To solve the system on a rectangular domain $[a_{x},b_{x}]\times[a_{y},b_{y}]$ we adopt uniform grid

$x_{i}=x_{0}+ih_{x},∼{}∼{}h_{x}=\frac{b_{x}-a_{x}}{N_{x}+1},∼{}∼{}y_{j}=y_{0}+% jh_{y},∼{}∼{}h_{y}=\frac{b_{y}-a_{y}}{N_{y}+1},∼{}∼{}i=0,\dots,N_{x},j=0,\dots% ,N_{y}.$

Figure 1.

Densities of Interacting Cells at $t=0.05,\alpha=1$ .

For simplicity, we denote explicit discrete Caputo fractional derivatives of $c_{1},c_{2},c_{3},u_{1},u_{2},v$ by ${\Delta}^{\alpha}_{\tau}c_{1},{\Delta}^{\alpha}_{\tau}c_{2},{\Delta}^{\alpha}_% {\tau}c_{3},{\Delta}^{\alpha}_{\tau}u_{1},{\Delta}^{\alpha}_{\tau}u_{2},{% \Delta}^{\alpha}_{\tau}v$ , i.e.

$\displaystyle{\Delta}^{\alpha}_{\tau}c^{n+1}_{1i,j}=\frac{1}{{\tau}^{\alpha}}% \left(c^{n+1}_{1i,j}-\sum^{k+1}_{n=1}{{\textrm{\^{a}}}^{\alpha}_{n}c^{k+1-n}_{% 1i,j}}\right)-{\tilde{c}}^{n}_{1i,j},$ (95) $\displaystyle{\Delta}^{\alpha}_{\tau}c^{n+1}_{2i,j}=\frac{1}{{\tau}^{\alpha}}% \left(c^{n+1}_{2i,j}-\sum^{k+1}_{n=1}{{\beta}^{\alpha}_{n}c^{k+1-n}_{2i,j}}% \right)-{\tilde{c}}^{n}_{2i,j},$ (96) $\displaystyle{\Delta}^{\alpha}_{\tau}c^{n+1}_{3i,j}=\frac{1}{{\tau}^{\alpha}}% \left(c^{n+1}_{3i,j}-\sum^{k+1}_{n=1}{{\beta}^{\alpha}_{n}c^{k+1-n}_{3i,j}}% \right)-{\tilde{c}}^{n}_{3i,j},$ (97) $\displaystyle{\Delta}^{\alpha}_{\tau}u^{n+1}_{1i,j}=\frac{1}{{\tau}^{\alpha}}% \left(u^{n+1}_{1i,j}-\sum^{k+1}_{n=1}{{\beta}^{\alpha}_{n}u^{k+1-n}_{1i,j}}% \right)-{\tilde{u}}^{n}_{1i,j},$ (98) $\displaystyle{\Delta}^{\alpha}_{\tau}u^{n+1}_{2i,j}=\frac{1}{{\tau}^{\alpha}}% \left(u^{n+1}_{2i,j}-\sum^{k+1}_{n=1}{{\beta}^{\alpha}_{n}u^{k+1-n}_{2i,j}}% \right)-{\tilde{u}}^{n}_{2i,j},$ (99) $\displaystyle{\Delta}^{\alpha}_{\tau}v^{n+1}_{i,j}=\frac{1}{{\tau}^{\alpha}}% \left(v^{n+1}_{i,j}-\sum^{k+1}_{n=1}{{\beta}^{\alpha}_{n}v^{k+1-n}_{i,j}}% \right)-{\tilde{v}}^{n}_{i,j}.$ (100)

For notational convenience, central in space discrete forms of 2 dimensions space derivatives is denoted by

$\displaystyle\Delta U^{n}_{i,j}=\frac{U^{n}_{i+1,j}-2U^{n}_{i,j}+U^{n}_{i-1,j}% }{h^{2}_{x}}+\frac{U^{n}_{i,j+1}-2U^{n}_{i,j}+U^{n}_{i,j-1}}{h^{2}_{y}},$ $\displaystyle\nabla U^{n}_{i,j}=\frac{U^{n}_{i+1,j}-U^{n}_{i-1,j}}{2h_{x}}+% \frac{U^{n}_{i,j+1}-U^{n}_{i,j-1}}{2h_{y}}.$ (101)

Figure 2.

Densities of Interacting Cells at $t=0.15,\alpha=1$ .

Now, substituting Eqs (95)–(5) into Eqs (23)–(28), we have the discretize system thus

$\displaystyle{\Delta}^{\alpha}_{\tau}c^{n+1}_{1i,j}=D_{1}\Delta c^{n}_{1i,j}+{% \rho}_{1}+b_{3}c^{n}_{3i,j}-{b_{1}c}^{n}_{1i,j}*v^{n}_{i,j}-{\gamma}_{1}c^{n}_% {1i,j},$ (102) $\displaystyle{\Delta}^{\alpha}_{\tau}c^{n+1}_{2i,j}=D_{2}\Delta c^{n}_{2i,j}-% \chi(c^{n}_{2i,j}*\Delta u^{n}_{1i,j}+\nabla c^{n}_{2i,j}*\nabla u^{n}_{1i,j})% +b_{1}c^{n}_{1i,j}*v^{n}_{i,j}{-}(b_{2}+{\gamma}_{2})c^{n}_{2i,j},$ (103) $\displaystyle{\Delta}^{\alpha}_{\tau}c^{n+1}_{3i,j}=D_{3}\Delta c^{n}_{3i,j}+b% _{2}c^{n}_{2i,j}-(b_{3}+{\gamma}_{3})c^{n}_{3i,j},$ (104) $\displaystyle{\Delta}^{\alpha}_{\tau}u^{n+1}_{1i,j}=D_{4}\Delta u^{n}_{1i,j}+{% \rho}_{2}-k_{3}c^{n}_{2i,j}*u^{n}_{1i,j}-k_{1}u^{n}_{1i,j}*v^{n}_{i,j}-{\mu}_{% 1}u^{n}_{1i,j},$ (105) $\displaystyle{\Delta}^{\alpha}_{\tau}u^{n+1}_{2i,j}=D_{5}\Delta u^{n}_{2i,j}+k% _{1}u^{n}_{1i,j}*v^{n}_{i,j}+k_{3}c^{n}_{2i,j}*u^{n}_{1i,j}-(k_{2}+{\mu}_{2})u% ^{n}_{2i,j},$ (106) $\displaystyle{\Delta}^{\alpha}_{\tau}v^{n+1}_{i,j}=D_{6}\Delta v^{n}_{i,j}+Nk_% {2}u^{n}_{2i,j}-{\mu}_{3}v^{n}_{i,j},$ (107)

where * is element wise multiplication of matrices. The stability of the scheme is ensured by

${\tau}^{\alpha}\leqslant\frac{h^{2}_{x}h^{2}_{y}}{2D_{6}(h^{2}_{x}+h^{2}_{y})}.$ (108)

Figure 3.

Density of Partially Mature DCs at $t=0.06,\alpha=0.5,\alpha=0.55,\alpha=0.6$ .

Figure 4.

Density of Infected T Cells at $t=0.06,\alpha=0.5,\alpha=0.55,\alpha=0.6$ .

6. Results and discussions

In this section, the results will be presented both for normal and subdiffusion cases. Recall that the normal diffusion system is recovered for $\alpha=1$ . In both cases, we used the parameter values in Table 1 and the following initial conditions

$\displaystyle c^{0}_{1}(x,y)=2000e^{-{0.5(x-1)}^{2}-{0.5(y-1)}^{2}},∼{}∼{}c^{0% }_{2}(x,y)=0,c^{0}_{3}(x,y)=0,$ $\displaystyle u^{0}_{1}(x,y)=3600e^{-{0.5(x+1)}^{2}-{0.5(y+1)}^{2}},∼{}∼{}u^{0% }_{2}(x,y)=0,v^{0}(x,y)=1.$ (109)

Using a square domain $[-4,4]\times[-4,4],$ it is assumed that the DCs have highest initial density at $(x,y)=(1,1)$ while susceptible T cells are densed at $(x,y)=(-1,-1).$ Notice that, before infection $c^{0}_{2}(x,y)=0,c^{0}_{3}(x,y)=0,u^{0}_{2}(x,y)=0,$ after which a viral particle was introduced, i.e. $v^{0}(x,y)=1$ . Furthermore, we take the time step size $h=0.0001$ and $N_{x}=N_{y}=100$ so that $h_{x}=h_{y}$ .

6.1 Normal diffusion

In Figs 1 and 2 above, we examine the cell densities at two different time scenarios for normal diffusion ( $\alpha=1$ ). In this case, a wave is observed in the plots of partially matured DCs (sub-figures 1b and 2b ) towards the T Cells which reflects also on the density of matured DCs. This implies increase in number of infected T cells with time increases as shown in Fig. 2e.

6.2 Subdiffusion

Here, the aim of the numerical study is to examine the impact of subdiffusion diffusion on the infection of T Cells. The plots are only shown for partially matured (activated) DCs (Fig. 3) and the density of infected T cells (Fig. 4). We inferred from the numerical experiments that, as the diffusion of DCs gets slower, the probability of the captured viral particles by DCs to infect T Cells becomes high. This reflects on density of infected T Cell shown in Fig. 4, as the time scaling parameter $\alpha\longrightarrow 1$ , the number of infected T cells decreases.

7. Conclusion

In this paper, a time-fractional diffusion model was proposed to study the effect of antigen presentation process by DCs on the progress of HIV. The spatial homogeneous equilibrium states of the model were obtain, the disease free case proved to be uniformly and asymptotically stable even with the presence of diffusion and chemotaxis. However, the endemic steady state is only stable with the presence of diffusion provided a threshold condition is fulfilled. A priori estimates were also obtained in the appropriate Sobolev spaces. The numerical simulations further revealed that, crowd induced slower diffusion (subdiffusion) affects the movement DCs and antigen presentation process in general. Unfortunately, subdiffusion due to crowded cells favours the pathogenesis of HIV by increasing the probability of infection of T cells during antigen presentation. Of course, this model can be generalized to include other immune cells and necessary factors. We recommend that, vivo experiments should investigate this result through cells tracking and other relevant technology. This may help the improvement of antiviral drugs and reduce chronic immune activation associated with the disease.

Footnotes

Acknowledgments

Authors are thankful to anonymous referees for their insightful inputs to increase the quality of manuscript.

References

Adams

and Fournier

, Sobolev spaces pure and applied mathematics (Amsterdam), 140 (2003).

Ammi

M.S.

and Taakili

, Finite difference method for the time-fractional thermistor problem, International Journal of Difference Equations 8(1) (2013), 77–97.

Arafa

Rida

and Khalil

, A fractional-order model of hiv infection with drug therapy effect, Journal of the Egyptian Mathematical Society 22(3) (2014), 538–543.

Bouchaud

J.-P.

and Georges

, Anomalous diffusion in disordered media: Statistical mechanisms, models and physical applications, Physics Reports 195(4-5) (1990), 127–293.

Brakken-Thal

, Gershgorin’s theorem for estimating eigenvalues, Source: http://buzzard.ups.edu/courses/2007spring/projects/brakkenthal-paper.pdf, 2007.

Brockmann

David

and Gallardo

A.M.

, Human mobility and spatial disease dynamics, Reviews of Nonlinear Dynamics and Complexity 2 (2009), 1–24.

Cartea

and del Castillo-Negrete

, Fractional diffusion models of option prices in markets with jumps, Physica A: Statistical Mechanics and Its Applications 374(2) (2007), 749–763.

Chun

T.-W.

Stuyver

Mizell

S.B.

Ehler

L.A.

, Mican

Baseler

Lloyd

A.L.

Nowak

M.A.

and Fauci

A.S.

, Presence of an inducible hiv-1 latent reservoir during highly active antiretroviral therapy, Proceedings of the National Academy of Sciences 94(24) (1997), 13193–13197.

Das

and Gupta

, A mathematical model on fractional lotka–volterra equations, Journal of Theoretical Biology 277(1) (2011), 1–6.

10.

del Castillo-Negrete

Carreras

and Lynch

, Nondiffusive transport in plasma turbulence: A fractional diffusion approach, Physical Review Letters 94(6) (2005), 065003.

11.

Demirci

Unal

and Özalp

, A fractional order seir model with density dependent death rate abstract full text, Hacettepe Journal of Mathematics and Statistics 40(2) (2011).

12.

Deng

Z.-Q.

De Lima

J.L.

De Lima

M.I.P.

and Singh

V.P.

, A fractional dispersion model for overland solute transport, Water Resources Research 42(3) (2006).

13.

Diethelm

, The analysis of fractional differential equations: An application-oriented exposition using differential operators of Caputo type, Springer, 2010.

14.

Diethelm

, A fractional calculus based model for the simulation of an outbreak of dengue fever, Nonlinear Dynamics 71(4) (2013), 613–619.

15.

Dimitrov

, Numerical approximations for fractional differential equations, arXiv preprint arXiv:1311.3935, 2013.

16.

Ding

and Ye

, A fractional-order differential equation model of hiv infection of cd4+ t-cells, Mathematical and Computer Modelling 50(3) (2009), 386–392.

17.

Flinsenberg

Spel

Jansen

Koning

de Haar

Plantinga

Scholman

van Loenen

Nierkens

Boon

et al., Cognate cd4 t-cell licensing of dendritic cells heralds anti-cytomegalovirus cd8 t-cell immunity after human allogeneic umbilical cord blood transplantation, Journal of Virology 89(2) (2015), 1058–1069.

18.

Granelli-Piperno

Delgado

Finkel

Paxton

and Steinman

R.M.

, Immature dendritic cells selectively replicate macrophagetropic (m-tropic) human immunodeficiency virus type 1, while mature cells efficiently transmit both m-and t-tropic virus to t cells, Journal of Virology 72(4) (1998), 2733–2737.

19.

Hanert

, Front dynamics in a two-species competition model driven by lévy flights, Journal of Theoretical Biology 300 (2012), 134–142.

20.

Hanert

Schumacher

and Deleersnijder

, Front dynamics in fractional-order epidemic models, Journal of Theoretical Biology 279(1) (2011), 9–16.

21.

Heath

S.L.

Tew

J.G.

Tew

J.G.

Szakal

A.K.

and Burton

G.F.

, Follicular dendritic cells and human immunodeficiency virus infectivity, Nature 377(6551) (1995), 740–744.

22.

D.D.

Neumann

A.U.

Perelson

A.S.

Chen

Leonard

J.M.

and Markowitz

, Rapid turnover of plasma virions and cd4 lymphocytes in hiv-1 infection, Nature 373(6510) (1995), 123–126.

23.

Höfling

and Franosch

, Anomalous transport in the crowded world of biological cells, Reports on Progress in Physics 76(4) (2013), 046602.

24.

Hogue

I.B.

Bajaria

S.H.

Fallert

B.A.

Qin

Reinhart

T.A.

and Kirschner

D.E.

, The dual role of dendritic cells in the immune response to human immunodeficiency virus type 1 infection, Journal of General Virology 89(9) (2008), 2228–2239.

25.

Kamath

A.T.

Pooley

O’Keeffe

M.A.

Vremec

Zhan

Lew

A.M.

D’Amico

Tough

D.F.

and Shortman

, The development, maturation, and turnover rate of mouse spleen dendritic cell populations, The Journal of Immunology 165(12) (2000), 6762–6770.

26.

Lanzavecchia

and Sallusto

, Regulation of t cell immunity by dendritic cells, Cell 106(3) (2001), 263–266.

27.

Lifson

J.D.

Rossio

J.L.

Piatak

Parks

Kiser

Coalter

Fisher

Flynn

B.M.

Czajak

et al., Role of cd8+ lymphocytes in control of simian immunodeficiency virus infection and resistance to rechallenge after transient early antiretroviral treatment, Journal of Virology 75(21) (2001), 10187–10199.

28.

Lizana

Ambjörnsson

Taloni

Barkai

and Lomholt

M.A.

, Foundation of fractional langevin equation: Harmonization of a many-body problem, Physical Review E 81(5) (2010), 051118.

29.

Lutz

Kollmann

and Bechinger

, Single-file diffusion of colloids in one-dimensional channels, Physical Review Letters 93(2) (2004), 026001.

30.

Mathai

A.M.

and Haubold

H.J.

, Special functions for applied scientists, Springer, 2008.

31.

Mempel

T.R.

Scimone

M.L.

Mora

J.R.

and von Andrian

U.H.

, In vivo imaging of leukocyte trafficking in blood vessels and tissues, Current Opinion in Immunology 16(4) (2004), 406–417.

32.

Metzler

and Klafter

, The random walk’s guide to anomalous diffusion: A fractional dynamics approach, Physics Reports 339(1) (2000), 1–77.

33.

Miller

M.J.

Wei

S.H.

Cahalan

M.D.

and Parker

, Autonomous t cell trafficking examined in vivo with intravital two-photon microscopy, Proceedings of the National Academy of Sciences 100(5) (2003), 2604–2609.

34.

Monel

Beaumont

Vendrame

Schwartz

Brand

and Mammano

, Hiv cell-to-cell transmission requires the production of infectious virus particles and does not proceed through env-mediated fusion pores, Journal of Virology 86(7) (2012), 3924–3933.

35.

Nowak

M.A.

Bangham

C.R.

et al., Population dynamics of immune responses to persistent viruses, Science-AAAS-Weekly Paper Edition 272(5258) (1996), 74–79.

36.

Odibat

Z.M.

and Momani

, An algorithm for the numerical solution of differential equations of fractional order, Journal of Applied Mathematics & Informatics 26(1_2) (2008), 15–27.

37.

Okada

Miller

M.J.

Parker

Krummel

M.F.

Neighbors

Hartley

S.B.

O’Garra

Cahalan

M.D.

and Cyster

J.G.

, Antigen-engaged b cells undergo chemotaxis toward the t zone and form motile conjugates with helper t cells, PLoS Biology 3(6) (2005), e150.

38.

Pawelek

K.A.

Liu

Pahlevani

and Rong

, A model of hiv-1 infection with two time delays: Mathematical analysis and comparison with patient data, Mathematical Biosciences 235(1) (2012), 98–109.

39.

Perelson

A.S.

Neumann

A.U.

Markowitz

Leonard

J.M.

D.D.

et al., Hiv-1 dynamics in vivo: Virion clearance rate, infected cell life-span, and viral generation time, Science 271(5255) (1996), 1582–1586.

40.

Petrás

, Fractional derivatives, fractional integrals, and fractional differential equations in matlab, In Engineering Education and Research Using MATLAB, InTech, 2011.

41.

Pinto

C.M.

and Carvalho

A.R.

, A latency fractional order model for hiv dynamics, Journal of Computational and Applied Mathematics 312 (2017), 240–256.

42.

Pletcher

R.H.

Tannehill

J.C.

and Anderson

, Computational fluid mechanics and heat transfer, CRC Press, 2012.

43.

Podlubny

, Fractional differential equations: An introduction to fractional derivatives, fractional differential equations, to methods of their solution and some of their applications, volume 198, Academic Press, 1998.

44.

Rosenbloom

D.I.

Hill

A.L.

Rabi

S.A.

Siliciano

R.F.

and Nowak

M.A.

, Antiretroviral dynamics determines hiv evolution and predicts therapy outcome, Nature Medicine 18(9) (2012), 1378–1385.

45.

Sardar

Rana

Bhattacharya

Al-Khaled

and Chattopadhyay

, A generic model for a single strain mosquito-transmitted disease with memory on the host and the vector, Mathematical Biosciences 263 (2015), 18–36.

46.

Scherer

Kalla

S.L.

Tang

and Huang

, The grünwald–letnikov method for fractional differential equations, Computers & Mathematics with Applications 62(3) (2011), 902–917.

47.

Smith

B.A.

Gartner

Liu

Perelson

A.S.

Stilianakis

N.I.

Keele

B.F.

Kerkering

T.M.

Ferreira-Gonzalez

Szakal

A.K.

Tew

J.G.

et al., Persistence of infectious hiv on follicular dendritic cells, The Journal of Immunology 166(1) (2001), 690–696.

48.

Stilianakis

N.I.

Boucher

De Jong

M.D.

Van Leeuwen

Schuurman

and De Boer

, Clinical data sets of human immunodeficiency virus type 1 reverse transcriptase-resistant mutants explained by a mathematical model, Journal of Virology 71(1) (1997), 161–168.

49.

Usero

, Fractional taylor series for caputo fractional derivatives. construction of numerical schemes, Dpto deMatemáticaAplicada, Universidad Complutense de Madrid, Spain, 2008.

50.

Wang

and Liu

, The stability analysis of a general viral infection model with distributed delays and multi-staged infected progression, Communications in Nonlinear Science and Numerical Simulation 20(1) (2015), 263–272.

51.

Wang

Zhang

and Kuniya

, Mathematical analysis for an age-structured hiv infection model with saturation infection rate, Electron J Differ Equ 33 (2015), 2015.

52.

Wei

Ghosh

S.K.

Taylor

M.E.

Johnson

V.A.

Emini

E.A.

Deutsch

Lifson

J.D.

Bonhoeffer

Nowak

M.A.

Hahn

B.H.

et al., Viral dynamics in human immunodeficiency virus type 1 infection, Nature 373(6510) (1995), 117–122.

53.

Wodarz

, Killer cell dynamics: Mathematical and computational approaches to immunology, 2006.

54.

Wodarz

Arnaout

R.A.

Nowak

M.A.

and Lifson

J.D.

, Transient antiretroviral treatment during acute simian immunodeficiency virus infection facilitates long–term control of the virus, Philosophical Transactions of the Royal Society of London B: Biological Sciences 355(1400) (2000), 1021–1029.

55.

and Zhou

, Bifurcation analysis of hiv-1 infection model with cell-to-cell transmission and immune response delay, Mathematical Biosciences and Engineering: MBE 13(2) (2016), 343–367.

56.

Yakimovich

Gumpert

Burckhardt

C.J.

Lütschg

V.A.

Jurgeit

Sbalzarini

I.F.

and Greber

U.F.

, Cell-free transmission of human adenovirus by passive mass transfer in cell culture simulated in a computer model, Journal of Virology 86(18) (2012), 10123–10137.

57.

Yang

Zou

and Ruan

, Global dynamics of a delayed within-host viral infection model with both virus-to-cell and cell-to-cell transmissions, Mathematical Biosciences 270 (2015), 183–191.

58.

Gao

and Ding

, A generalized gronwall inequality and its application to a fractional differential equation, Journal of Mathematical Analysis and Applications 328(2) (2007), 1075–1081.

59.

Zhang

Huang

Liu

and Fan

, Stability analysis for a fractional hiv infection model with nonlinear incidence, Discrete Dynamics in Nature and Society (2015), 2015.

60.

Zhuang

and Liu

, Implicit difference approximation for the time fractional diffusion equation, Journal of Applied Mathematics and Computing 22(3) (2006), 87–99.