
Citation: Emmanuel Reyes-Uribe, Nathalia Serna-Marquez, Eduardo Perez Salazar. DDRs: receptors that mediate adhesion, migration and invasion in breast cancer cells[J]. AIMS Biophysics, 2015, 2(3): 303-317. doi: 10.3934/biophy.2015.3.303
[1] | Jinhu Xu, Yicang Zhou . Bifurcation analysis of HIV-1 infection model with cell-to-cell transmission and immune response delay. Mathematical Biosciences and Engineering, 2016, 13(2): 343-367. doi: 10.3934/mbe.2015006 |
[2] | Yu Yang, Gang Huang, Yueping Dong . Stability and Hopf bifurcation of an HIV infection model with two time delays. Mathematical Biosciences and Engineering, 2023, 20(2): 1938-1959. doi: 10.3934/mbe.2023089 |
[3] | Shunyi Li . Hopf bifurcation, stability switches and chaos in a prey-predator system with three stage structure and two time delays. Mathematical Biosciences and Engineering, 2019, 16(6): 6934-6961. doi: 10.3934/mbe.2019348 |
[4] | Dan Liu, Shigui Ruan, Deming Zhu . Stable periodic oscillations in a two-stage cancer model of tumor and immune system interactions. Mathematical Biosciences and Engineering, 2012, 9(2): 347-368. doi: 10.3934/mbe.2012.9.347 |
[5] | Jiawei Deng, Ping Jiang, Hongying Shu . Viral infection dynamics with mitosis, intracellular delays and immune response. Mathematical Biosciences and Engineering, 2023, 20(2): 2937-2963. doi: 10.3934/mbe.2023139 |
[6] | Kaushik Dehingia, Anusmita Das, Evren Hincal, Kamyar Hosseini, Sayed M. El Din . Within-host delay differential model for SARS-CoV-2 kinetics with saturated antiviral responses. Mathematical Biosciences and Engineering, 2023, 20(11): 20025-20049. doi: 10.3934/mbe.2023887 |
[7] | Yan Wang, Minmin Lu, Daqing Jiang . Viral dynamics of a latent HIV infection model with Beddington-DeAngelis incidence function, B-cell immune response and multiple delays. Mathematical Biosciences and Engineering, 2021, 18(1): 274-299. doi: 10.3934/mbe.2021014 |
[8] | Wanxiao Xu, Ping Jiang, Hongying Shu, Shanshan Tong . Modeling the fear effect in the predator-prey dynamics with an age structure in the predators. Mathematical Biosciences and Engineering, 2023, 20(7): 12625-12648. doi: 10.3934/mbe.2023562 |
[9] | Huan Kong, Guohong Zhang, Kaifa Wang . Stability and Hopf bifurcation in a virus model with self-proliferation and delayed activation of immune cells. Mathematical Biosciences and Engineering, 2020, 17(5): 4384-4405. doi: 10.3934/mbe.2020242 |
[10] | Juan Wang, Chunyang Qin, Yuming Chen, Xia Wang . Hopf bifurcation in a CTL-inclusive HIV-1 infection model with two time delays. Mathematical Biosciences and Engineering, 2019, 16(4): 2587-2612. doi: 10.3934/mbe.2019130 |
According to clinical data, drug treatment against persistent human infections sometimes fails to consistently eradicate the infection from the host [1,2,3], such as human immunodeficiency virus (HIV), hepatitis B virus (HBV), and hepatitis C virus (HCV). For all these viruses, drug treatment is not effective, and for HIV, lifelong therapy is commonly required to control viral replication [4,5]. The main reason is that these viruses are able to suppress and impair immune responses, resulting in persistent infection [6,7]. An alternative strategy is to use drug therapies to boost virus-specific immune response and then induce sustained viral suppression in the absence of drugs. Antiviral therapy aimed at boosting specific immune responses has attracted more and more attentions recently. In 2003, Komarova et al. [6] used mathematical models to study immune response dynamics during therapy in the context of immunosuppressive infections. They showed that a single phase of antiviral drug therapy is also possible to establish sustained immunity. They assumed that the degree of immune expansion depends on virus load and the response inhibits virus growth. This assumption is natural for any branch of the adaptive immune system, such as CD4 T cells, CD8 T cells, or antibodies. Denote by $ y(t) $ and $ z(t) $, respectively, the virus population and immune cell population at time $ t $. Komarova et al. [6] proposed the following system of ordinary differential equations:
$ y′(t)=ry(t)(1−y(t)K)−py(t)z(t)−by(t),z′(t)=cy(t)z(t)1+dy(t)−qy(t)z(t)−μz(t). $ | (1.1) |
The virus population grows at a rate described by the logistic function. The viral replication rate is a linear decreasing function of viral loads with a maximum value of $ r $ and it vanishes at a viral load $ K $. The treatment is modeled as a reduction in $ r $. Immune cells kill virus at a rate $ pyz $. $ b $ and $ \mu $ are the natural decay rates of viral particles and immune cells, respectively. Immune expansion is modeled by the growth function $ {cyz\over 1+dy} $. When the virus load is low, the level of immune response is simply proportional to both the viral load and the immune cells. The immune response saturates when the virus load is sufficiently high. Immune cells are inhibited by the virus at a rate $ qyz $. Komarova et al. presented a simple relationship between the timing of therapy and efficacy of the drugs. In the presence of strong viral suppression, they showed that therapy should be stopped relatively early because a longer duration of treatment may lead to a failure. On the other hand, in the presence of weaker viral suppression, stopping treatment too early is detrimental, and therapy has to be continued beyond a time threshold. Essentially, model (1.1) has two stable equilibria: a virus dominant equilibrium with no sustained immunity and an immune control equilibrium with sustained immunity. This bistability allows a solution from the basin of the attraction of the virus dominant equilibrium to be lifted to that of the immune control equilibrium via a single phase of therapy.
Model (1.1) assumes that the immune response is instantaneous. Nevertheless, it has been established that the immune response process involves a sequence of events such as antigenic activation, selection, and proliferation of the immune cells [8]. Moreover, oscillatory viral loads and immune cells were observed from clinical data [9], while the oscillatory behavior is not exhibited in model (1.1). A natural question is what is the cause of the oscillations. In 2014, Shu et al. [10] incorporated the time lag during the immune response process into model (1.1). By studying the dynamics between an immunosuppressive infection and antiviral immune response, they demonstrated that the time lag is the main factor causing oscillations. Here, $ \tau $ is denoted as the time lag for the immune system to trigger a sequence of events. The activation rate of immune cells at time $ t $ depends on the virus load and the number of immune cells at time $ t-\tau $. The above assumptions lead to the following system:
$ y′(t)=ry(t)(1−y(t)K)−py(t)z(t)−by(t),z′(t)=cy(t−τ)z(t−τ)1+dy(t−τ)−qy(t)z(t)−μz(t). $ | (1.2) |
It was shown that the delayed antiviral immune response may induce sustained periodic oscillations, transient oscillations and even sustained aperiodic oscillations (chaos), and the model can admit two stable equilibria and can also allow a stable equilibrium to coexist with a stable periodic solution.
In the aforementioned work, the mortality of the population during the time lag has been ignored. Consideration of the survival probability during the time lag requires an additional delay-dependent multiplier in the nonlinear delayed feedback term. The death of the free viral particles depends on the mature immune cells instead of all the newly generated immune cells. We denote $ z(t) $ as the population of mature immune cells at time $ t $, and consider the following delay differential equation with a delay-dependent coefficient:
$ y′(t)=f(y(t))−py(t)z(t)−by(t),z′(t)=e−δτh(y(t−τ))z(t−τ)−qy(t)z(t)−μz(t). $ | (1.3) |
Here, $ f(y) $ describes the virus population growth function; $ \delta > 0 $ is the death rate of the immature immune cells during the time lag; $ e^{-\delta\tau} $ represents the survival probability of the immature immune cells surviving $ \tau $ time units before becoming mature; $ h(y)z $ denotes the immune cell growth function; the other parameters have the same meaning as in model (1.2). In addition to more general virus and immune cell growth functions, model (1.3) differs from model (1.2) in the sense that it introduces an additional delay-dependent survival probability $ e^{-\delta\tau} $. As we shall see later, this quantity will make a significant difference in the bifurcation analysis. For instance, model (1.3) does not possess any positive equilibrium for sufficiently large $ \tau $ and its global Hopf branches always have bounded $ \tau $ components. However, for model (1.2), the existence condition of positive equilibrium is independent of $ \tau $ and the global Hopf branches always have unbounded $ \tau $ components. Furthermore, there exists only one stability switch of positive equilibrium for model (1.2), but the stability of positive equilibrium for our model (1.3) may switch multiple times.
Actually, the model (1.3) can be derived from a stage structured population model for $ u(t, a) $ as below
$ y′(t)=f(y(t))−py(t)z(t)−by(t),∂tu(t,a)+∂au(t,a)=−(μ(a)+q(a)y(t))u(t,a), $ |
with the stage-specific decay rate of immune cells, $ \mu(a) $, inhibition rate of immune cells, $ q(a) $,
$ \mu(a) = \left\{μ,a≥τ,δ,a<τ,\right. \ \ \ \ q(a) = \left\{q,a≥τ,0,a<τ,\right. $ |
where $ u(t, a) $ is the population of immune cells at age $ a $ and time $ t $. Note that the population of mature immune cells is $ z(t) = \int_{\tau}^{\infty}u(t, a) da $. We integrate along characteristic lines to find
$ z′(t)=u(t,τ)−u(t,∞)−qy(t)z(t)−μz(t)=u(t−τ,0)e−δτ−u(t,∞)−qy(t)z(t)−μz(t). $ |
Note that $ u(t, \infty) = 0 $ and $ u(t-\tau, 0) = h(y(t-\tau))z(t-\tau) $. The equation for $ z(t) $ follows.
Throughout this paper, we assume that the virus population growth function $ f(y) $ and the immune expansion rate $ h(y) $ satisfy the following conditions:
(H1) $ f(y) $ is continuously differentiable, $ f(0) = 0 $ and $ f'(0) > b $; there exists $ \bar{y} > 0 $ such that $ f(\bar{y}) = 0 $, and $ (y-\bar{y})f(y) < 0 $ for positive $ y\neq\bar{y} $; $ f''(y) < 0 $ for $ y\in[0, \bar{y}] $.
(H2) $ h(y) $ is continuously differentiable, $ h(0) = 0 $ and $ h'(0) > q $; $ h'(y) > 0 $ and $ h''(y) < 0 $ for all $ y\geq0 $.
The assumptions (H1)-(H2) are biological conditions for infection dynamics. For instance, $ f(0) = 0 $ means no new infection in the absence of virus; $ f'(0) > b $ indicates that the intrinsic growth rate is greater than the decay rate. Similarly, $ h(0) = 0 $ implies that immune cells cannot reproduce without the presence of virus; $ h'(0) > q $ guarantees that the initial expansion rate is greater than the inhibition rate.
The rest of this paper is organized as follows. We first present some preliminary results concerning the well-posedness of model (1.3), and describe global dynamics of the ordinary differential equation system (2.4) in Section 2. Local and global stability analysis of equilibria as well as local and global Hopf bifurcation analysis are given in Section 3. In Section 4, numerical simulations based on the bifurcation analysis are reported and discussed. Finally, we give a brief summary in Section 5.
To establish the well-posedness of model (1.3), we choose the phase space $ \mathcal{C}\times\mathcal{C} $, where $ \mathcal{C} $ is the Banach space of continuous functions on $ [-\tau, 0] $ defined by $ \mathcal{C}: = \{\phi:[-\tau, 0]\rightarrow \mathbb{R} \ \text{is continuous} \} $, and the norm is defined by $ \|\phi\| = \sup\limits_{-\tau\leq\theta\leq0}\phi(\theta) $. The nonnegative cone of $ \mathcal{C} $ is denoted by $ \mathcal{C^{+}} $. As usual, $ \phi_t\in\mathcal{C} $ is defined by $ \phi_t(\theta) = \phi(t+\theta) $ for $ \theta\in[-\tau, 0] $. For biological applications, the initial condition of (1.3) is given as
$ (y_0,z_0)\in X: = \mathcal{C}^{+}\times\mathcal{C}^{+}. $ |
The existence and uniqueness of the solution of model (1.3) follows from the standard theory of functional differential equations [11]. Using the same manner as in proving [12,Proposition 2.1], we can easily show that the solution of (1.3) with initial conditions in $ X $ is nonnegative, which implies that $ X $ is positively invariant under the solution map of (1.3).
Next, we shall prove that all solutions of (1.3) are ultimately bounded. It follows from the first equation in (1.3) that $ y'(t)\leq f(y) $, which yields $ \limsup\limits_{t\rightarrow\infty}y(t)\leq\bar y $. By the nonnegativity of the solution of (1.3) and $ (\textbf{H}_{1}) $-$ (\textbf{H}_{2}) $, we obtain
$ (y(t)+ph′(0)z(t+τ))′≤Mf+(e−δτh(y)h′(0)y−1)py(t)z(t)−by(t)−μph′(0)z(t+τ)≤Mf−γ(y(t)+ph′(0)z(t+τ)), $ |
where $ M_f = \max\limits_{y\geq0}f(y) $ and $ \gamma = \min\{b, \mu\} $. Thus,
$ \limsup\limits_{t\rightarrow\infty}\left(y(t)+\frac{p}{h'(0)}z(t+\tau)\right)\leq {M_f\over \gamma}, $ |
which implies that $ y(t) $ and $ z(t) $ are ultimately bounded in $ X $. To summarize, we obtain the following proposition.
Proposition 2.1. The solutions of model (1.3) with the initial conditions in $ X $ are nonnegative, and the region
$ \Gamma = \left\{(\phi, \psi)\in X: \|\phi\|\leq \bar{y}, \ \phi(-\tau)+{p\over h'(0)}\psi(0)\leq {M_f\over \gamma}\right\} $ |
is positively invariant and absorbing in $ X $; namely, all solutions ultimately enter $ \Gamma $.
Clearly, model (1.3) always has a trivial equilibrium $ E_0 = (0, 0) $. Since $ f'(0) > b $, there also exists a boundary equilibrium $ E_b = (\tilde{y}, 0) $, where $ \tilde{y} $ is the unique positive root of $ f(y) = by $. We call $ E_b $ the virus dominant equilibrium (VDE). Assume that $ E^* = (y^*, z^*) $ is a positive equilibrium with $ y^* > 0 $ and $ z^* > 0 $, then it must satisfy the following equilibrium equations
$ e−δτh(y∗)=qy∗+μ, z∗=f(y∗)−by∗py∗. $ | (2.1) |
Obviously, $ z^* > 0 $ if and only if $ y^* < \tilde{y} $. Hence, $ E^* $ exists if and only if $ y^* $ is a positive root of
$ h(y)y+μ/q=qeδτ $ | (2.2) |
with $ y^* < \tilde{y} $. It is easily to calculate that $ \text{sign}\left({h(y)\over y+\mu/q}\right)' = \text{sign} \ g(y) $, where $ g(y) = h'(y)(y+{\mu\over q})-h(y) $. Note that $ g(0) > 0, g(\infty) < 0 $ and $ g'(y) = h''(y)(y+\mu/q) < 0 $. There exists a unique $ y_c > 0 $ such that $ g(y_c) = 0 $, and $ (y-y_c)g(y) < 0 $ for $ y\neq y_c $. Therefore, $ {h(y)\over y+\mu/q} $ is strictly increasing on $ [0, y_c) $, and strictly decreasing on $ (y_c, \infty) $. We denote
$ s:=h(yc)yc+μ/q=supy≥0h(y)y+μ/q. $ | (2.3) |
If $ s < qe^{\delta\tau} $, (2.2) has no positive root. If $ s\ge qe^{\delta\tau} $, (2.2) has exactly two positive roots (counting multiplicity) $ y_1^* $ and $ y_2^* $ such that $ y_1^*\le y_c\le y_2^* $. These two roots coincide (i.e., $ y_1^* = y_2^* = y_c $) if and only if $ s = qe^{\delta\tau} $. Throughout this paper, we make the following assumption to ensure the existence of positive roots for (2.2).
(H3) $ s > q $ and $ \tau\le\tau_{max} $, where $ \tau_{max} = \frac{1}{\delta}\ln\frac s q $.
If ($ \textbf{H}_{3} $) is violated, there does not exist any positive equilibrium for our model and the boundary equilibrium is locally asymptotically stable. Summarizing the above analysis, we obtain the following result.
Proposition 2.2. Assume that (H1)-(H3) are satisfied.
(i) If $ y_1^*\geq \tilde{y} $ holds, then there are two equilibria: $ E_{0} = (0, 0) $ and $ E_{b} = (\tilde{y}, 0) $.
(ii) If $ y_1^* < \tilde{y}\leq y_2^* \ \mathit{\text{or}} \ \ y_1^* = y_2^* < \tilde{y} $ holds, then there are three equilibria: $ E_{0}, $ $ E_{b} $ and $ E_1^* = (y_1^*, z_1^*) $.
(iii) If $ y_1^* < y_2^* < \tilde{y} $ holds, then there are four equilibria: $ E_{0}, $ $ E_b, $ $ E_1^* $ and $ E_2^* = (y_2^*, z_2^*). $
When $ \tau = 0 $, model (1.3) reduces to an ODE system which generalizes (1.1)
$ y′(t)=f(y(t))−py(t)z(t)−by(t),z′(t)=h(y(t))z(t)−qy(t)z(t)−μz(t). $ | (2.4) |
We next provide a complete description about its global dynamics.
Lemma 2.3. System (2.4) has no closed orbits.
Proof. Using an argument similar to Proposition 2.1, we can show that
$ \Gamma_0 = \left\{(\phi, \psi)\in \mathbb{R}_+^2: \phi\leq \bar{y}, \ \phi+{p\over h'(0)}\psi\leq {M_f\over \gamma}\right\} $ |
is positively invariant and absorbing in $ \mathbb{R}_+^2 $. Thus, it suffices to show that System (2.4) has no closed orbits in $ \Gamma_0 $.
Let $ (P(y, z), Q(y, z)) $ denote the vector field of (2.4), then we have
$ \frac{\partial(BP)}{\partial y}+\frac{\partial(BQ)}{\partial z} = \frac{1}{z(t)}(\frac{f(y)}{y})', \ \ \text{where} \ \ B(y, z) = {1\over yz}. $ |
Note that ($ \textbf{H}_{1} $) implied that $ (\frac{f(y)}{y})' < 0 $ for all $ y\leq \bar y $. Thus, $ \frac{\partial(BP)}{\partial y}+\frac{\partial(BQ)}{\partial z} < 0 $ in $ \Gamma_0 $. The nonexistence of closed orbits for (2.4) follows from the classical Dulac-Bendixson criterion. This ends the proof.
We linearize system (2.4) at $ E_0 = (0, 0) $ and calculate that the two eigenvalues are $ f'(0)-b > 0 $ and $ -\mu < 0 $. Thus, $ E_0 $ is a saddle point. For the linearized system at $ E_b = (\tilde{y}, 0) $, the two eigenvalues are
$ \lambda_1 = f'(\tilde{y})-b \lt 0 \ \ \text{and} \ \ \lambda_2 = h(\tilde{y})-q\tilde{y}-\mu. $ |
From the proof of Proposition 2.2, we know $ h(\tilde{y})-q\tilde{y}-\mu < 0 $ if $ \tilde{y} < y_1^* $ or $ \tilde{y} > y_2^* $, and $ h(\tilde{y})-q\tilde{y}-\mu > 0 $ provided that $ y_1^* < \tilde{y} < y_2^* $. This shows that $ E_b $ is stable if $ y_1^* > \tilde{y} $ or $ y_2^* < \tilde{y} $, and $ E_b $ is a saddle point if $ y_1^* < \tilde{y} < y_2^* $. For the linearized system at $ E_1^* = (y_1^*, z_1^*) $, the associated characteristic equation is
$ \lambda^2+(-f'(y_1^*)+pz_1^*+b)\lambda+py_1^*z_1^*(h'(y_1^*)-q) = 0. $ |
It follows from (2.1) and ($ \textbf{H}_{1} $) that $ -f'(y_1^*)+pz_1^*+b = -f'(y_1^*)+f(y_1^*)/y_1^* > 0. $ This, together with $ h'(y_1^*)-q > 0 $, implies that all eigenvalues must have negative real parts, and hence $ E_1^* $ is asymptotically stable. At $ E_2^* = (y_2^*, z_2^*) $, the two eigenvalues $ \lambda_1 $ and $ \lambda_2 $ satisfy
$ \lambda_1+\lambda_2 = f'(y_2^*)-{f(y_2^*)\over y_2^*} \lt 0, \ \ \lambda_1\lambda_2 = py_2^*z_2^*(h'(y_2^*)-q) \lt 0. $ |
This implies that one eigenvalue must be positive and the other negative. Thus, $ E_2^* $ is a saddle point whenever it exists. Summarizing the above analysis, we obtain the following global stability results on model (2.4).
Theorem 2.4. Assume that ($ \textbf{H}_{1} $)-($ \textbf{H}_{2} $) are satisfied and $ s > q $.
(i) If $ y_1^* > \tilde{y} $ holds, then $ E_0 $ is a saddle point and $ E_b $ attracts all solutions with initial conditions in $ \{(y(t), z(t))\in \mathbb{R}_+^2: y(0) > 0, z(0)\geq0\} $.
(ii) If $ y_1^* < \tilde{y} < y_2^* $ holds, then both $ E_0 $ and $ E_b $ are saddle points and $ E_1^* $ attracts all solutions with initial conditions in $ \{(y(t), z(t))\in \mathbb{R}_+^2: y(0) > 0, z(0) > 0\} $.
(iii) If $ y_1^* < y_2^* < \tilde{y} $ holds, then both $ E_0 $ and $ E_2^* $ are saddle points, and for any initial condition $ (y(0), z(0)) $ with $ y(0) > 0 $, the solution of (2.4) approaches either $ E_b $ or $ E_1^* $, the stable manifold of $ E_2^* $ separate the basins of attraction of $ E_b $ and $ E_1^* $.
In this section, we study the dynamics of model (1.3) with $ \tau > 0 $.
By linearizing (1.3) at $ E_0 = (0, 0) $, we find two eigenvalues $ f'(0)-b > 0 $ and $ -\mu < 0 $. Hence, $ E_0 $ is a saddle point. The characteristic equation associated with the linearization of model (1.3) at $ E_b $ is
$ (\lambda-f'(\tilde{y})+b)\left(\lambda+q\tilde{y}+\mu-e^{-\delta\tau}h(\tilde{y})e^{-\lambda\tau}\right) = 0. $ |
One eigenvalue is $ \lambda_{1} = f'(\tilde{y})-b < 0 $. Hence $ E_b $ is locally asymptotically stable if all zeros of
$ \Delta_b(\lambda) = \lambda+q\tilde{y}+\mu-e^{-\delta\tau}h(\tilde{y})e^{-\lambda\tau} $ |
have negative real parts. By [13,Lemma 6], there exists at least one positive eigenvalue if $ q\tilde{y}+\mu < e^{-\delta\tau}h(\tilde{y}) $, or equivalently, $ y_1^* < \tilde{y} < y_2^* $. On the other hand, if $ q\tilde{y}+\mu > e^{-\delta\tau}h(\tilde{y}) $, that is, either $ \tilde{y} < y_1^* $ or $ \tilde{y} > y_2^* $, then all eigenvalues have negative real parts. For the critical case $ q\tilde{y}+\mu = e^{-\delta\tau}h(\tilde{y}) $; namely, either $ \tilde{y} = y_1^* $ or $ \tilde{y} = y_2^* $, one eigenvalue is $ 0 $, and all other eigenvalues have negative real parts. By using the method of normal forms, we obtain that $ E_b $ is locally asymptotically stable. To summarize, we have the following result on the local stability of $ E_0 $ and $ E_b $.
Theorem 3.1. Consider model (1.3) under the assumptions ($ \textbf{H}_{1} $)-($ \textbf{H}_{3} $). $ E_0 $ is a saddle point. If either $ \tilde{y}\le y_1^* $ or $ \tilde{y}\ge y_2^* $, then $ E_b $ is locally asymptotically stable; while $ E_b $ is unstable if $ y_1^* < \tilde{y} < y_2^* $.
To establish the global stability of $ E_b $, we first show that the infection is persistent.
Lemma 3.2. $ \liminf\limits_{t\rightarrow\infty}y(t) > 0 $ for any solution of (1.3) with the initial condition in $ X_1 = \{(\phi, \psi)\in X: \phi(0) > 0, \psi(0)\geq0\} $.
Proof. We claim that $ \limsup\limits_{t\rightarrow\infty}y(t) > 0 $. If not, then $ y(t)\to0 $ as $ t\to\infty $. It then follows from (1.3) that $ z(t)\to0 $ as $ t\to\infty $. This contradicts to the fact that $ E_0 $ is unstable. Thus, $ y(t) $ is weakly persistent. By using [14,Theorem 3.4.6] and Proposition 2.1, there exists a global attractor for the solution semiflow of (1.3). This, together with [15,Theorem 2.3], implies that $ y(t) $ is actually strongly persistent; that is, $ \liminf\limits_{t\rightarrow\infty}y(t) > 0 $.
Our next result is concerned with global stability of $ E_b $.
Theorem 3.3. Assume that ($ \textbf{H}_{1} $)-($ \textbf{H}_{3} $) hold. Then $ E_b $ of model (1.3) is globally asymptotically stable in the region $ X_1 $ provided that $ \tau\geq\tau_b $, where $ \tau_b: = {1\over \delta}\ln{h'(0)\tilde{y}\over q\tilde{y}+\mu} $. Moreover, the condition $ \tau\geq\tau_b $ implies that $ \tilde{y}\leq y_1^* $.
Proof. By Proposition 2.1 and Theorem 3.1, it suffices to show that $ E_b $ is globally attractive in $ X_1\cap\Gamma $, which is a positively invariant set of (1.3). Let $ (y(t), z(t)) $ be a solution of (1.3) with the initial condition in $ X_1 $, it follows from Lemma 3.2 that $ y(t) > 0 $ for $ t > 0 $. Motivated by the earlier work in [16], we construct a Lyapunov functional $ L: X_1\cap\Gamma\rightarrow\mathbb{R} $ by
$ L(y_t,z_t) = y_t(0)-\tilde{y}\ln(y_t(0))+cz_t(0)+ce^{-\delta\tau}\int^{0}_{-\tau}h(y_{t}(\theta))z_{t}(\theta)d\theta, $ |
where $ c $ is a constant to be determined later. Calculating the time derivative of $ L $ along the solution, and using $ b = f(\tilde{y})/\tilde{y} $, $ h(y)\leq h'(0)y $ for all $ y\geq0 $, we obtain
$ dLdt=f(y)y(y−˜y)−b(y−˜y)+(p˜y−cμ)z+ce−δτh(y)z−(cq+p)yz,≤(f(y)y−f(˜y)˜y)(y−˜y)+(p˜y−cμ)z+(ce−δτh′(0)−cq−p)yz. $ |
Assumption ($ \textbf{H}_{1} $) implies that $ (\frac{f(y)}{y}-\frac{f(\tilde{y})}{\tilde{y}})(y-\tilde{y})\leq0 $ for all $ y\leq\bar y $. Since $ h'(0)e^{-\delta\tau}-q > 0 $, the condition $ \tau\geq\tau_b $ is the same as $ {p\tilde y\over\mu}\leq {p\over e^{-\delta\tau}h'(0)-q} $. We may choose a constant $ c\in[{p\tilde y\over\mu}, {p\over e^{-\delta\tau}h'(0)-q}] $ such that $ dL/dt\leq0 $. Moreover, $ dL/dt = 0 $ if and only if $ y(t) = \tilde y $ and $ z(t) = 0 $. Thus the maximal compact invariant set in $ \{dL/dt = 0\} $ is the singleton $ \{E_b\} $. By the LaSalle invariance principle [11], $ E_b $ is globally attractive in $ X_1\cap\Gamma $. Since $ X_1\cap\Gamma $ is absorbing in $ X_1 $, we conclude that $ E_b $ is globally attractive in $ X_1 $. By Theorem 3.2, $ E_b $ is globally asymptotically stable in $ X_1 $ if $ \tau\geq\tau_b $.
If $ \tau\geq\tau_b $, then $ e^{-\delta\tau}h'(0)\tilde y-q\tilde y-\mu\leq0 = e^{-\delta\tau}h(y_1^*)-qy_1^*-\mu \le e^{-\delta\tau}h'(0)y_1^*-qy_1^*-\mu $. Since $ e^{-\delta\tau}h'(0) > q $, we obtain $ \tilde y\le y_1^* $.
Here, we remark that if $ y_1^* = y_2^* < \tilde{y} $, then two positive equilibria $ E_1^* $ and $ E_2^* $ coincide, and $ E_b $ is locally asymptotically stable.
In this subsection, we investigate the stability of $ E_1^* $ and identify parameter regions in which the time delay can destabilize $ E_1^* $, lead to Hopf bifurcation and induce sustained oscillations.
The characteristic equation associated with the linearization of system (1.3) at equilibrium $ E_i^* \ (i = 1, 2) $ is
$ Δi(λ)=λ2+a1,iλ+a0,i+(b1,iλ+b0,i)e−λτ=0, $ | (3.1) |
where
$ a0,i=(f(y∗i)y∗i−f′(y∗i))h(y∗i)e−δτ−pqy∗iz∗i, a1,i=f(y∗i)y∗i−f′(y∗i)+h(y∗i)e−δτ>0,b0,i=(f′(y∗i)−f(y∗i)y∗i)h(y∗i)e−δτ+py∗iz∗ih′(y∗i)e−δτ, b1,i=−h(y∗i)e−δτ<0. $ |
We first investigate the stability of $ E_1^* $, regarding $ \tau $ as the bifurcation parameter. In the proof of Theorem 2.4, we have demonstrated that, when $ \tau = 0 $, all eigenvalues of (3.1) with $ i = 1 $ lie to the left of the imaginary axis and $ E_1^* $ is locally asymptotically stable. As $ \tau $ increases, $ E_1^* $ may lose its stability only when some eigenvalues cross the imaginary axis to the right. In view of
$ a_{0,1}+b_{0,1} = p y_{1}^{*}z_{1}^{*}(h'(y_{1}^{*})e^{-\delta\tau}-q) \gt 0, $ |
$ 0 $ is not an eigenvalue for any $ \tau\geq0 $. For simplicity, we drop the subscript $ 1 $ in the following arguments. Substituting $ \lambda = i\omega (\omega > 0) $ into (3.1) and separating the real and imaginary parts, we obtain
$ ω2−a0=b0cos(ωτ)+b1ωsin(ωτ),−a1ω=b1ωcos(ωτ)−b0sin(ωτ). $ |
Squaring and adding both the above equations lead to
$ G(ω,τ):=ω4+c1(τ)ω2+c0(τ)=0, $ | (3.2) |
where
$ c_1(\tau) = ({f(y_{1}^{*})\over y_{1}^{*}}-f'(y_{1}^{*}))^{2}+2pqy_{1}^{*}z_{1}^{*} \gt 0, \ \ c_0(\tau) = a_{0}^{2}-b_{0}^{2}. $ |
Then $ i\omega \; (\omega > 0) $ is an imaginary root of (3.1) if and only if $ G(\omega, \tau) = 0 $. Since $ c_1(\tau) > 0 $ for all $ 0\leq\tau\le\tau_{max} $. we know that $ G(\omega, \tau) = 0 $ has exactly one positive root if and only if $ c_0(\tau) < 0 $, or equivalently, $ a_0 < b_0 $.
If $ a_0 < b_0 $, the implicit function theorem implies that there exists a unique $ C^1 $ function $ \omega = \omega^*(\tau) $ such that $ G(\omega^*(\tau), \tau) = 0 $ for $ 0 < \tau\le\tau_{max} $. For $ i\omega^*(\tau) $ to be a root of (3.1), $ \omega^*(\tau) $ needs to satisfy the system
$ sin(ω(τ)τ)=b1ω(τ)3+(a1b0−a0b1)ω(τ)b21ω(τ)2+b20:=g1(τ),cos(ω(τ)τ)=(b0−a1b1)ω(τ)2−a0b0b21ω(τ)2+b20:=g2(τ). $ | (3.3) |
Set
$ I={τ:0≤τ≤τmax satisfies a0<b0}. $ | (3.4) |
For $ \tau\in I $, let $ \theta(\tau) $ be the unique solution of $ \sin\theta = g_1 $ and $ \cos\theta = g_2 $ in $ (0, 2\pi] $; that is,
$ \theta(\tau) = \left\{ arccosg2(τ), if ω(τ)2≤a0−a1b0b1,2π−arccosg2(τ), if ω(τ)2>a0−a1b0b1. \right. $ |
Following Beretta and Kuang [17], we define
$ Sn(τ)=τ−θ(τ)+2nπω(τ) for τ∈I with n∈N. $ | (3.5) |
One can check that $ \pm i\omega^*(\tau^*) $ are a pair of purely imaginary eigenvalues of (3.1) if and only if $ \tau^* $ is a zero of function $ S_n(\tau) $ for some $ n\in\mathbb{N} $. From [17,Theorem 2.2], we have
$ \text{Sign}\left({d Re\lambda(\tau)\over d\tau}\big|_{\tau = \tau^*}\right) = \text{Sign}\left({\partial G\over\partial\omega}(\omega^*(\tau^*),\tau^*)\right) \text{Sign}\; S'_n(\tau^*). $ |
Note that $ {\partial G\over\partial\omega}(\omega^*(\tau^*), \tau^*) > 0 $, thus we have the following result concerning the transversality condition.
Lemma 3.4. Assume that $ S_n(\tau) $ has a positive root $ \tau^*\in I $ for some $ n\in\mathbb{N} $, then a pair of simple purely imaginary roots $ \pm i\omega^*(\tau^*) $ of (3.1) exist at $ \tau = \tau^* $, and
$ \mathit{\text{Sign}}\left({d Re\lambda(\tau)\over d\tau}\big|_{\tau = \tau^*}\right) = \mathit{\text{Sign}}\; S'_n(\tau^*). $ |
Moreover, this pair of simple purely imaginary roots $ \pm i\omega^*(\tau^*) $ cross the imaginary axis from left to right at $ \tau = \tau^* $ if $ S'_n(\tau^*) > 0 $, and from right to left if $ S'_n(\tau^*) < 0 $.
It is easily seen that $ S_n(\tau) > S_{n+1}(\tau) $ for all $ \tau\in I, n\in\mathbb{N} $. When $ \tau = 0 $, the asymptotic stability of $ E_1^* $ implies that $ S_0(0) < 0 $. Denote
$ ˆτ=supI:=sup{τ:0≤τ≤τmax satisfies a0<b0}. $ | (3.6) |
If $ \hat\tau < \tau_{max} $, then $ c_0(\hat\tau) = 0 $, and thus $ \omega(\tau)\rightarrow0 $ as $ \tau\rightarrow\hat\tau^-. $ This, together with (3.3), yields $ \lim\limits_{\tau\rightarrow\hat\tau^-}\sin\theta(\tau) = 0 $ and $ \lim\limits_{\tau\rightarrow\hat\tau^-}\cos\theta(\tau) = -1 $. Therefore, $ \lim\limits_{\tau\rightarrow\hat\tau^-}\theta(\tau) = \pi, \ \lim\limits_{\tau\rightarrow\hat\tau^-}S_n(\tau) = -\infty. $
If $ \sup\limits_{\tau\in I} S_0(\tau) < 0 $, $ S_n(\tau) $ has no zeros in $ I $ for all $ n\in\mathbb{N} $, which excludes the existence of purely imaginary eigenvalues and thus implies that $ E_1^* $ is locally asymptotically stable for all $ \tau\in[0, \tau_{max}] $.
If $ \sup\limits_{\tau\in I} S_0(\tau) = 0 $, $ S_0(\tau) $ has a double zero in $ I $, denoted by $ \tau_* $, and $ S_0'(\tau_*) = 0 $. This, together with Lemma 3.4, implies that the transversality condition at $ \tau_* $ is not satisfied and all eigenvalues remain to the left of the pure imaginary axis. Thus, $ E_1^* $ is still locally asymptotically stable for $ \tau\in[0, \tau_{max}] $.
If $ \sup\limits_{\tau\in I} S_0(\tau) > 0 $, it then follows from $ S_0(0) < 0 $ and $ \lim\limits_{\tau\rightarrow\hat\tau^-}S_n(\tau) = -\infty $ that $ S_0(\tau) $ has at least two zeros in $ I $. For simplicity, we assume that
(H4) $ \sup\limits_{\tau\in I} S_0(\tau) > 0 $ and all zeros of $ S_n(\tau) $ with odd multiplicity are simple.
Under assumption ($ \textbf{H}_{4} $), each $ S_n(\tau) $ has even number of simple zeros. Now, we collect all simple zeros of $ S_n(\tau) $ with $ n\ge0 $ and list them in increasing order: $ 0 < \tau_0 < \tau_1 < \cdots < \tau_{2K-1} < \hat\tau $ $ (K\in\mathbb{N}^+) $. For each integer $ 0\leq i\leq K-1 $, we apply Lemma 3.4 to obtain $ S'_m(\tau_i) > 0 $ and $ S'_m(\tau_{2K-i-1}) < 0 $ for some $ m\in\mathbb{N} $. Therefore, the pair of simple conjugate purely imaginary eigenvalues $ \pm i\omega(\tau_i) $ crosses the imaginary axis from left to right, and the pair of simple conjugate purely imaginary eigenvalues $ \pm i\omega(\tau_{2K-i-1}) $ crosses the imaginary axis from right to left. Thus, system (1.3) undergoes a Hopf bifurcation at $ E_1^* $ when $ \tau = \tau_j \ (0\le j\le 2K-1) $. Moreover, $ E_1^* $ is asymptotically stable for $ \tau\in[0, \tau_0)\cup(\tau_{2K-1}, \tau_{max}] $, and unstable for $ \tau\in(\tau_0, \tau_{2K-1}) $.
For each $ k = 0, \cdots, 2K-1 $, there exists $ n_k $ such that $ \tau_k $ is a simple zero of $ S_{n_k}(\tau) $. Let $ T_k $ be the period of periodic solution bifurcated at $ \tau_k $. Applying the Hopf bifurcation theorem for delay differential equations [11,18], we have
$ T_k = {2\pi\over \omega(\tau_k)} = {2\pi\tau_k\over \theta(\tau_k)+2n_k\pi}. $ |
This, together with $ \theta(\tau_k)\in(0, 2\pi] $, implies that $ T_k > \tau_k $ if $ n_k = 0 $, and $ {\tau_k\over n_k+1}\le T_k < {\tau_k\over n_k} $ if $ n_k > 0 $. To summarize, we have the following results on stability of $ E_1^* $ and Hopf bifurcation.
Theorem 3.5. Assume that ($ \textbf{H}_{1} $)-($ \textbf{H}_{3} $) hold, $ y_1^* < \tilde y $ and $ \hat\tau < \tau_{max} $. Let $ I, S_n(\tau), \tau_{max}, \hat\tau $ be defined in (3.4), (3.5), ($ \textbf{H}_{3} $) and (3.6), respectively.
(i) If either $ I = \emptyset $ or $ \sup\limits_{\tau\in I} S_0(\tau)\le 0 $, then $ E_1^* $ is locally asymptotically stable for all $ \tau\in[0, \tau_{max}] $.
(ii) If ($ \textbf{H}_{4} $) holds, then there exist exactly $ 2K $ local Hopf bifurcation values, namely, $ 0 < \tau_0 < \tau_1 < \cdots < \tau_{2K-1} < \hat\tau $ such that model (3.1) undergoes a Hopf bifurcation at $ E_1^* $ when $ \tau = \tau_j $ for $ 0\le j\le 2K-1 $. $ E_1^* $ is locally asymptotically stable for $ \tau\in[0, \tau_0)\cup(\tau_{2K-1}, \tau_{max}] $, and unstable for $ \tau\in(\tau_0, \tau_{2K-1}) $. Moreover, the period $ T_k $ of periodic solution bifurcated at $ \tau_k $ satisfies $ T_k > \tau_k $ if $ \tau_k $ is a simple zero of $ S_0(\tau) $, and $ {\tau_k\over n_k+1}\le T_k < {\tau_k\over n_k} $ if $ \tau_k $ is a simple zero of $ S_{n_k}(\tau) $ with $ n_k > 0 $.
Theorem 3.6. Assume that ($ \textbf{H}_{1} $)-($ \textbf{H}_{3} $) are satisfied. If $ y_1^* < y_2^* < \tilde{y} $ holds, then $ E_2^* $ is unstable for all $ \tau\geq 0 $. Moreover, if $ a_{0, 2}\le b_{0, 2} $, then the characteristic equation (3.1) with $ i = 2 $ has no purely imaginary eigenvalues.
Proof. Note that in (3.1) with $ i = 2 $, $ \Delta_2(0) = p y_{2}^{*}z_{2}^{*}(h'(y_{2}^{*})e^{-\delta\tau}-q) < 0 $. For any $ \tau\geq 0 $, we know that $ \Delta_2(\infty) > 0 $, thus as long as $ E_2^* $ exists, the associated characteristic equation must have at least one positive real root and hence $ E_2^* $ is always unstable for $ \tau\geq 0 $.
Using a similar argument in the study of root distribution for (3.1) at $ E_1^* $, we note that $ c_{1, 2}(\tau) > 0, c_{0, 2}(\tau) = a_{0, 2}^2-b_{0, 2}^2 $ and $ a_{0, 2}+b_{0, 2} = p y_{2}^{*}z_{2}^{*}(h'(y_{2}^{*})e^{-\delta\tau}-q) < 0 $. Thus, $ G(\omega, \tau) $ in (3.2) has no positive zeros if $ a_{0, 2}\le b_{0, 2} $. Therefore, (3.1) with $ i = 2 $ has no purely imaginary eigenvalues if $ a_{0, 2}\le b_{0, 2} $.
Since $ E_2^* $ is unstable (and thus biologically irrelevant) whenever it exists, we name $ E_1^* $, instead of $ E_2^* $, the immune control equilibrium (ICE).
Note that Theorem 3.5 states that if $ y_1^* < \tilde y $ and ($ \textbf{H}_{4} $) holds, then periodic solutions can bifurcate from $ E_1^* $ when $ \tau $ is near the local Hopf bifurcation values $ \tau_k, k = 0, 1, 2, \ldots, 2K-1. $ For integer $ n\ge0 $, we denote
$ J:={τ0,τ1,⋯,τ2K−1}, Jn:={τ∈J:Sn(τ)=0}, J0=J∖J0. $ | (3.7) |
According to Theorem 3.5, $ J_n $ has even number of bifurcation values, and the period of any periodic solution bifurcated near a bifurcation point in $ J^0 $ is bounded from both below and above, while the periodic solution bifurcated near a bifurcation point in $ J_0 $ has a lower bound for its period. As we shall see later in the numerical simulation, it seems impossible to find an upper bound for such period. In what follows, we will restrict our investigation on the set $ J^0 $ and assume that $ J^0\neq\emptyset $. Especially, we will discuss the global continuation of periodic solutions bifurcated from the point $ (E_1^*, \tau_k) $ with $ \tau_k\in J^0 $ as the bifurcation parameter $ \tau $ varies. We shall use the global Hopf bifurcation theorem for delay differential equations [19] and show that model (1.3) admits periodic solutions globally for all $ \tau\in(\underline\tau, \bar\tau) $, where $ \underline\tau = \min J^0 $ and $ \bar\tau = \max J^0 $.
Let $ x(t) = (y(\tau t), z(\tau t))^T $, model (1.3) can be rewritten as a general functional differential equation:
$ x′(t)=F(xt,τ,T),(t,τ,T)∈R+×I×R+, $ | (3.8) |
where $ x_{t}(\theta) = x(t+\theta) $ for $ \theta\in[-1, 0] $, and $ x_{t}\in \mathcal{X}: = C([-1, 0], \mathbb{R}_{+}^2) $, and
$ F(xt,τ,T)=(τf(x1t(0))−pτx1t(0)x2t(0)−bτx1t(0)τe−δτh(x1t(−1))x2t(−1)−qτx1t(0)x2t(0)−μτx2t(0)) $ | (3.9) |
with $ x_t = (x_{1t}, x_{2t})\in\mathcal{X} $. Identifying the subspace of $ \mathcal{X} $ consisting of all constant functions with $ \mathbb{R}_{+}^2 $, we obtain a restricted function
$ ˜F(x,τ,T):=F∣R2+×I×R+=(τf(x1)−pτx1x2−bτx1τe−δτh(x1)x2−qτx1x2−μτx2). $ |
Obviously, $ \widetilde{F} $ is twice continuously differentiable, i.e., the condition $ (\textbf{A1}) $ in [19] holds. We denote the set of stationary solutions of system (3.8) by $ {\bf{E}}(F) = \{(\widetilde{x}, \widetilde{\tau}, \widetilde{T})\in\mathbb{R}_{+}^2\times I\times\mathbb{R}_{+}: \; \widetilde{F}(\widetilde{x}, \widetilde{\tau}, \widetilde{T}) = 0 \}. $ It follows from Proposition 2.2 that (ⅰ) if $ y_1^* < \tilde{y} < y_2^* $ holds,
$ {\bf{E}}(F) = \{(E_0,\tau,T), (E_b,\tau,T), (E_1^*,\tau,T); \; (\tau,T)\in I\times\mathbb{R}_{+}\}; $ |
(ⅱ) if $ y_1^* < y_2^* < \tilde{y} $ holds,
$ {\bf{E}}(F) = \{(E_0,\tau,T), (E_b,\tau,T), (E_1^*,\tau,T), (E_2^*,\tau,T); \; (\tau,T)\in I\times\mathbb{R}_{+}\}. $ |
For any stationary solution $ (\widetilde{x}, \tau, T) $, the characteristic function is
$ \Delta_{(\widetilde{x},\tau,T)}(\lambda) = \lambda \text{Id}-DF(\widetilde{x},\tau,T)(e^{\lambda\cdot}\text{Id}). $ |
Note that if either $ y_1^* < \tilde{y} < y_2^* $ or $ y_1^* < y_2^* < \tilde{y} $ holds, $ 0 $ is not a eigenvalue of any stationary solution of (3.8) and hence the condition (A2) in [19] holds. It can be checked easily from (3.9) that the smoothness condition (A3) in [19] is satisfied.
Theorem 3.5 implies that if $ y_1^* < \tilde y $ and ($ \textbf{H}_{4} $) holds, then for each integer $ 0\le k\le 2K-1 $ the stationary solution $ (E_1^*, \tau_k, 2\pi/(\omega_k\tau_k)) $ is an isolated center of (3.8), where $ \omega_k = \omega(\tau_k) $ is the unique positive zero of $ G(\omega, \tau) $ in (3.2), and there is only one purely imaginary characteristic value of the form $ im (2\pi/\widetilde{T}) $ with $ m = 1 $ and $ \widetilde{T} = 2\pi/(\omega_k\tau_k) $. Moreover, if follows from Lemma 3.4 that the crossing number $ \zeta_{1}(E_1^*, \tau_k, 2\pi/(\omega_k\tau_k)) $ at each of these centers is
$ ζ1(E∗1,τk,2πωkτk)=−Sign(dReλ(τ)dτ|τ=τk)={−1, 0≤k≤K−1,1, K≤k≤2K−1. $ | (3.10) |
Thus the condition (A4) in [19] holds. We then define a closed subset $ \Sigma(F) $ of $ \mathcal{X}\times I \times \mathbb{R}_+ $ by
$ \Sigma(F) = \mathcal{C}l\{(x,\tau,T)\in \mathcal{X}\times I \times \mathbb{R}_+: x \mbox{ is a nontrivial T-periodic soltuion of (3.8)}\}, $ |
and for each integer $ 0\le k\le 2K-1 $, let $ {\bf{C}}(E_1^*, \tau_k, 2\pi/(\omega_k\tau_k)) $ denote the connected component of $ (E_1^*, \tau_k, 2\pi/(\omega_k\tau_k)) $ in $ \Sigma(F) $. By Theorem 3.5, $ {\bf{C}}(E_1^*, \tau_k, 2\pi/(\omega_k\tau_k)) $ is a nonempty subset of $ \Sigma(F) $.
To find the interval of $ \tau $ in which periodic solutions exist, we shall further investigate the properties of periodic solutions of (3.8).
Lemma 3.7. All nonconstant and nonnegative periodic solutions of (3.8) are uniformly bounded. Actually, we have $ 0 < y(t), z(t)\le M $ for any $ t\in\mathbb{R} $, where $ M = \max\{\bar y, h'(0)M_f/(p\gamma)\} $.
Proof. Since $ y(t) $ is nonconstant and nonnegative, there exists $ t_0 > 0 $ such that $ y(t_0) > 0 $. Integrating the equation for $ y'(t) $ gives
$ e^{\int_{t_0}^t[pz(s)+b]ds}y(t) = y(t_0)+\int_{t_0}^t e^{\int_{t_0}^r[pz(s)+b]ds}f(y(r))dr \gt 0,\; t \gt t_0. $ |
Hence, $ y(t) > 0 $ for all $ t > t_0 $. Since $ y $ is periodic, we have $ y(t) > 0 $ for all $ t > 0 $. Similarly, if $ z(t_0) > 0 $ for some $ t_0 > 0 $, we integrate the equation for $ z'(t) $ to obtain
$ e^{\int_{t_0}^t[qy(s)+\mu]ds}z(t) = z(t_0)+\int_{t_0}^t e^{\int_{t_0}^r[qy(s)+\mu]ds}e^{-\delta\tau}h(y(r-\tau))z(r-\tau)dr \gt 0 $ |
for all $ t > t_0 $. This together with periodicity of $ z $ implies that $ z(t) > 0 $ for all $ t > 0 $.
By Proposition 2.1, we have $ \limsup\limits_{t\rightarrow\infty}y(t)\le\bar y $. We claim $ y(t)\le M $ for all $ t\ge 0 $. Otherwise, if $ y(t_1) > M $ for some $ t_1\ge0 $, then
$ \lim\limits_{n\rightarrow\infty}y(t_1+nT) = y(t_1) \gt M, $ |
where $ T $ is the period of the nonconstant and nonnegative periodic solutions. This contradicts the fact that $ y(t) $ is eventually bounded above by $ M $ as $ t\rightarrow\infty $. Thus, $ M $ is the upper bound of $ y(t) $. Similarly, from Proposition 2.1, we have $ \limsup\limits_{t\rightarrow\infty}z(t)\le h'(0)M_f/(p\gamma)\le M $. This implies that $ z(t) $ is uniformly bounded by $ M $.
Lemma 3.8. System (3.8) has no nontrivial periodic solution of period 1.
Proof. Note that any nontrivial $ 1 $-periodic solution $ x(t) $ of system (3.8) is also a nontrivial periodic solution of model (2.4), which does not admit any nontrivial periodic solution via Lemma 2.3.
We are now in the position to state our result concerning the properties of the global Hopf branches.
Theorem 3.9. Assume that ($ \textbf{H}_{1} $)-($ \textbf{H}_{4} $), either $ y_1^* < \tilde{y} < y_2^* $ or $ y_1^* < y_2^* < \tilde{y} $ holds, $ \hat\tau < \tau_{max}, J^0\neq\emptyset $ and $ a_{0, 2}\le b_{0, 2} $. Then for system (3.8), we have the following results:
(i) $ $The connected component $ {\bf{C}}(E_1^*, \tau_k, 2\pi/(\omega_k\tau_k)) $ is bounded for all $ \tau_k\in J^0 $.
(ii) Let $ n\ge1 $. If $ J_n $ has $ 2k_n $ bifurcation values, ordered as $ \tau_{n, 1} < \cdots < \tau_{n, 2k_n} $, then for each $ i = 1, \cdots, k_n $, there exists $ j > k_n $ such that $ \tau_{n, i} $ and $ \tau_{n, j} $ are connected by a global Hopf branch. Similarly, for each $ j = k_n+1, \cdots, 2k_n $, there exists $ i\le k_n $ such that $ \tau_{n, i} $ and $ \tau_{n, j} $ are connected by a global Hopf branch. Especially, if $ k_n = 1 $, then the two bifurcation values of $ J_n $ are connected.
(iii) For each $ \tau\in(\min J_n, \max J_n) $ with integer $ n\ge1 $, system (3.8) has at least one periodic solution with period in $ (1/(n+1), 1/n) $.
Proof. Lemma 3.7 implies that the projection of $ {\bf{C}}(E_1^*, \tau_k, 2\pi/(\omega_k\tau_k)) $ onto $ \mathcal{X} $ is bounded. Note that system (3.8) has no periodic solutions of period $ 1 $, and thus no periodic solutions of period $ 1/n_{k} $ or $ 1/(n_{k}+1) $ for any positive integer $ n_{k} $. It follows from Lemma 3.8 that the period $ T $ of a periodic solution on the connected component $ {\bf{C}}(E_1^*, \tau_k, 2\pi/(\omega_k\tau_k)) $ satisfies
$ {1\over n_{k}+1} \lt T \lt {1\over n_{k}} $ |
with $ \tau_{k}\in J_{n_k} $ and $ n_k\ge1 $. Hence, the projection of $ {\bf{C}}(E_1^*, \tau_k, 2\pi/(\omega_k\tau_k)) $ onto the $ T $-space is bounded for $ \tau_k\in J^0 $. Note that $ \tau\in I $ and $ I $ is a bounded interval. Therefore, $ {\bf{C}}(E_1^*, \tau_k, 2\pi/(\omega_k\tau_k)) $ is bounded in $ \mathcal{X}\times I\times \mathbb{R}_+ $ for any $ \tau_k\in J^0 $. This proves (ⅰ).
From the proof of Theorems 3.1 and 3.6, we know that the stationary solutions $ (E_0, \tau, T) $, $ (E_b, \tau, T) $ and $ (E_2^*, \tau, T) $ cannot be a center for any $ \tau $ and $ T $ if either $ y_1^* < \tilde{y} < y_2^* $ or $ y_1^* < y_2^* < \tilde{y} $ holds and $ a_{0, 2}\le b_{0, 2} $. It then follows from the global bifurcation theorem [19,Theorem 3.3] that $ \mathcal{E}: = {\bf{C}}(E_1^*, \tau_k, 2\pi/(\omega_k\tau_k))\cap {\bf{E}}(F) $ is finite and
$ \sum\limits_{(\widetilde{x},\tau,T)\in \mathcal{E}}\zeta_1(\widetilde{x},\tau,T) = 0. $ |
If $ J_n $ has $ 2k_n $ bifurcation values, ordered as $ \tau_{n, 1} < \cdots < \tau_{n, 2k_n} $, then by (3.10), $ \zeta_1(\widetilde{x}, \tau_{n, i}, T) = -1 $ for each $ i = 1, \cdots, k_n $, and $ \zeta_1(\widetilde{x}, \tau_{n, j}, T) = 1 $ for each $ j = k_n+1, \cdots, 2k_n $. This together with the above equation implies (ⅱ).
Note that $ \tau_{n, 1} = \min J_n $ and $ \tau_{n, 2k_n} = \max J_n $. The $ \tau $ projection of the global bifurcation branch bifurcated from $ \tau_{n, 1} $ includes $ (\tau_{n, 1}, \tau_{n, k_n+1}) $, while the $ \tau $ projection of the global bifurcation branch bifurcated from $ \tau_{n, 2k_n} $ includes $ (\tau_{n, k_n}, \tau_{n, 2k_n}) $. Thus, for any $ \tau\in(\tau_{n, 1}, \tau_{n, 2k}) $, system (3.8) has at least one periodic solution. According to the proof of (ⅰ), the period of this periodic solution lies in $ (1/(n+1), 1/n) $. Thus, (ⅲ) is proved.
In this section, we present some numerical simulations to demonstrate our theoretical results. We will use the Matlab package DDE-BIFTOOL developed by Engelborghs et al. [20,21] to sketch the global Hopf branches $ {\bf{C}}(E_1^*, \tau_k, 2\pi/(\omega_k\tau_k)) $ for $ 0\le k\le 2K-1 $. Following the study of ODE model in [6], we choose
$ f(y) = ry(1-{y\over K}), \ \ h(y) = {cy\over 1+dy}, $ |
and set the parameter values as follows:
$ r=6day−1,K=3virusmm−3,p=1mm3 cells−1 day−1,b=3day−1,c=4mm3 virus−1day−1,d=0.5mm3 virus−1,δ=0.03day−1,q=1mm3 virus−1day−1,μ=0.5day−1. $ | (4.1) |
Here, the units for viral loads $ y $ and immune cell density $ z $ are virus per $ \text{mm}^3 $ and cells per $ \text{mm}^3 $, respectively. The delay $ \tau $ has the unit in days. It is easy to calculate
$ \hat\tau = 17.1939 \lt \tau_{max} = 19.1788 \lt \tau_b = 36.6204. $ |
From Proposition 2.2, there are exactly two equilibria $ E_{0} = (0, 0) $ and $ E_{b} = (1.5, 0) $ when $ \tau > \tau_{max} $, and, if $ \tau\leq\tau_{max} $, there exists a third equilibrium $ E_1^* = (y_1^*, z_1^*) $, where $ y_1^*, z_1^* $ depend on $ \tau $. This is the unique positive equilibrium when $ \tau\in[0, \tau_e] $ and there exists another positive equilibrium $ E_2^* = (y_2^*, z_2^*) $ when $ \tau\in(\tau_e, \tau_{max}) $; see Figure 1. Here, $ \tau_e $ is the critical value when $ y_2^* = \tilde y $; that is,
$ \tau_e = {1\over\delta}\ln{c\tilde y\over(q\tilde y+\mu)(1+d\tilde y)} = 17.9666. $ |
At $ \tau = \tau_{max} $, the two positive equilibria coincide: $ y_1^* = y_2^* $.
Theorem 3.1 implies that $ E_0 $ is a saddle point, and $ E_b $ is locally asymptotically stable for $ \tau\ge\tau_e $, and unstable for $ \tau < \tau_e $, Moreover, by Theorem 3.3, $ E_b $ is globally asymptotically stable in the region $ X_1 $ provided that $ \tau\geq\tau_b $. It can be verified that the conditions of case (ⅱ) in Theorem 3.5 are satisfied only if $ 0\leq\tau < \hat\tau $. By Theorem 3.5, we know that there are exactly $ 4 $ local Hopf bifurcation values with $ K = 2 $, namely,
$ \tau_0\approx0.2994 \lt \tau_1\approx8.9274 \lt \tau_2\approx13.9287 \lt \tau_3\approx17.0322, $ |
as shown in Figure 2. Correspondingly,
$ \omega_0\approx0.9537 \gt \omega_1\approx0.7540 \gt \omega_2\approx0.4994 \gt \omega_3\approx0.1069. $ |
Clearly, $ E_1^* $ is stable for $ \tau\in[0, \tau_0)\cup(\tau_3, \tau_{max}] $ and unstable for $ \tau\in(\tau_0, \tau_3) $, which implies that $ \tau_0 $ and $ \tau_3 $ are stability switches. Remark that $ \tau_3 < \tau_e $.
In Figure 3, we plot the four global Hopf branches $ {\bf{C}}(E_1^*, \tau_k, 2\pi/(\omega_k\tau_k)) $ with $ 0\le k\le 3 $. It is observed that the branches $ {\bf{C}}(E_1^*, \tau_1, 2\pi/(\omega_1\tau_1)) $ and $ {\bf{C}}(E_1^*, \tau_2, 2\pi/(\omega_2\tau_2)) $ are bounded and connected, which agrees with Theorem 3.9. The periodic solutions on these global Hopf branches are plotted on the phase plane of $ y $ and $ z $; see Figures 4–6. As we know, the stability of the periodic solution is completely determined by the associated principal Floquet multiplier: if the principal Floquet multiplier is larger than $ 1 $, then the corresponding periodic solution is unstable, otherwise, the bifurcated periodic solution is stable [11].
By using DDE-BIFTOOL, we can further calculate the associated principal Floquet multipliers; see Figure 7. Note that the periodic solution on the first branch $ {\bf{C}}(E_1^*, \tau_0, 2\pi/(\omega_0\tau_0)) $ is stable for small $ \tau $, and becomes unstable as $ \tau $ increases; the periodic solution on the second branch $ {\bf{C}}(E_1^*, \tau_1, 2\pi/(\omega_1\tau_1)) $–or equivalently, the third branch $ {\bf{C}}(E_1^*, \tau_2, 2\pi/(\omega_2\tau_2)) $–is always unstable; the periodic solution on the fourth branch $ {\bf{C}}(E_1^*, \tau_3, 2\pi/(\omega_3\tau_3)) $ is stable for $ \tau $ near $ \tau_3 $, and becomes unstable as $ \tau $ decreases. Moreover, for any $ \tau\in(\tau_1, \tau_2) $, model (1.3) has at least one periodic solution on the second (or third) global Hopf branch with period in the interval $ (\tau/2, \tau) $, as shown in Figure 8.
It is interesting to note from Figure 3 that the first and the fourth Hopf branches do not connect. This phenomenon is quite different from that in the existing literature. We have proved that both the periodic solutions and delays on any global Hopf branch are uniformly bounded. According to the global Hopf bifurcation theory [19], either the first and the fourth Hopf branches connect, or these two branches have unbounded periods. From Figure 8, we observe that the periods of periodic oscillations on the fourth branch $ {\bf{C}}(E_1^*, \tau_3, 2\pi/(\omega_3\tau_3)) $ seem to be unbounded. To further explore this interesting phenomenon, we numerically sketch in Figure 9 a bifurcation diagram for model (1.3). It is suggested that model (1.3) can have chaotic solutions and thus allow aperiodic oscillations when $ \tau $ lies in the interval where the first and fourth Hopf branches disappear. Figure 9 also confirms that $ \tau_0 $ and $ \tau_3 $ are stability switches: (ⅰ) when $ \tau $ crosses $ \tau_0 $ from left to right, the positive equilibrium $ E_1^* $ becomes unstable, while a stable periodic solution bifurcates; and (ⅱ) when $ \tau $ crosses $ \tau_3 $ from left to right, $ E_1^* $ regains its stability while there exists a stable periodic solution when $ \tau $ is close to $ \tau_3 $ from the left. It should be mentioned that the bifurcation diagram has a jump at $ \tau_{max} $, which is due to the fact that there exists no positive equilibrium when $ \tau > \tau_{max} $ but the positive equilibrium at $ \tau = \tau_{max} $ is nonzero.
Considering the fact the immune response process involves a sequence of events such as antigenic activation, selection, and proliferation of the immune cells and the mortality of immune cells during the time lag, we derive an immunosuppressive infection model from a stage structured population model. We find that the time lag $ \tau $ plays a key role in the infection and immune response dynamics. As $ \tau $ increases, the model dynamics shifts through four possible outcomes: (ⅰ) virus persists without immune response when $ \tau > \tau_{max} $; (ⅱ) the outcome is initial condition dependent and delay dependent, either the immunity can be completely destroyed, or the immune mediated control can be established when $ \tau\in(\tau_e, \tau_{max}] $; (ⅲ) immune response controls the virus when $ \tau\in[0, \tau_0)\cup(\tau_{2K-1}, \tau_e) $; and (ⅳ) sustained oscillation appears or chaotic dynamics occurs when $ \tau\in(\tau_0, \tau_{2K-1}) $. Our theoretical results coincide with the phenomenon of oscillatory viral loads and immune cells observed in clinical data [9].
From numerical simulation, we also observe an unbounded global Hopf bifurcation branch with bounded periodic solutions, bounded delays, but unbounded periods. To the best of our knowledge, this phenomenon is new, because, in the existing literature, either the global Hopf branch is bounded or it is unbounded with unbounded delays. A detailed theoretical study of this interesting numerical phenomenon should enrich the global Hopf bifurcation theory for delay differential equations, and we leave it as a future work.
We are very grateful to the anonymous referees for their valuable comments and suggestions, which help to improve the presentation of this manuscript. H. Shu was partially supported by the National Natural Science Foundation of China (No.11971285, No.11601392), and the Fundamental Research Funds for the Central Universities.
The authors declare there is no conflict of interest.
[1] |
Gjorevski N, Nelson CM (2011) Integrated morphodynamic signalling of the mammary gland. Nat Rev Mol Cell Biol 12: 581-593. doi: 10.1038/nrm3168
![]() |
[2] |
Forsyth IA, Neville MC (2009) Introduction: hormonal regulation of mammary development and milk protein gene expression at the whole animal and molecular levels. J Mammary Gland Biol Neoplasia 14: 317-319. doi: 10.1007/s10911-009-9146-4
![]() |
[3] | Polyak K, Kalluri R (2010) The role of the microenvironment in mammary gland development and cancer. Cold Spring Harb Perspect Biol 2: a003244. |
[4] |
Vogel WF, Abdulhussein R, Ford CE (2006) Sensing extracellular matrix: an update on discoidin domain receptor function. Cell Signal 18: 1108-1116. doi: 10.1016/j.cellsig.2006.02.012
![]() |
[5] |
Ozbek S, Balasubramanian PG, Chiquet-Ehrismann R, et al. (2010) The evolution of extracellular matrix. Mol Biol Cell 21: 4300-4305. doi: 10.1091/mbc.E10-03-0251
![]() |
[6] |
Egeblad M, Rasch MG, Weaver VM (2010) Dynamic interplay between the collagen scaffold and tumor evolution. Curr Opin Cell Biol 22: 697-706. doi: 10.1016/j.ceb.2010.08.015
![]() |
[7] |
McNally S, Martin F (2011) Molecular regulators of pubertal mammary gland development. Ann Med 43: 212-234. doi: 10.3109/07853890.2011.554425
![]() |
[8] |
Hinck L, Silberstein GB (2005) Key stages in mammary gland development: the mammary end bud as a motile organ. Breast Cancer Res 7: 245-251. doi: 10.1186/bcr1331
![]() |
[9] | Watson CJ (2006) Post-lactational mammary gland regression: molecular basis and implications for breast cancer. Expert Rev Mol Med 8: 1-15. |
[10] |
Brisken C, Kaur S, Chavarria TE, et al. (1999) Prolactin controls mammary gland development via direct and indirect mechanisms. Dev Biol 210: 96-106. doi: 10.1006/dbio.1999.9271
![]() |
[11] |
Oakes SR, Rogers RL, Naylor MJ, et al. (2008) Prolactin regulation of mammary gland development. J Mammary Gland Biol Neoplasia 13: 13-28. doi: 10.1007/s10911-008-9069-5
![]() |
[12] |
Sternlicht MD, Kouros-Mehr H, Lu P, et al. (2006) Hormonal and local control of mammary branching morphogenesis. Differentiation 74: 365-381. doi: 10.1111/j.1432-0436.2006.00105.x
![]() |
[13] | Fata JE, Werb Z, Bissell MJ (2004) Regulation of mammary gland branching morphogenesis by the extracellular matrix and its remodeling enzymes. Breast Cancer Res 6: 1-11. |
[14] |
Myllyharju J, Kivirikko KI (2004) Collagens, modifying enzymes and their mutations in humans, flies and worms. Trends Genet 20: 33-43. doi: 10.1016/j.tig.2003.11.004
![]() |
[15] |
Yeh YC, Lin HH, Tang MJ (2012) A tale of two collagen receptors, integrin beta1 and discoidin domain receptor 1, in epithelial cell differentiation. Am J Physiol Cell Physiol 303: C1207-1217. doi: 10.1152/ajpcell.00253.2012
![]() |
[16] |
Carafoli F, Hohenester E (2013) Collagen recognition and transmembrane signalling by discoidin domain receptors. Biochim Biophys Acta 1834: 2187-2194. doi: 10.1016/j.bbapap.2012.10.014
![]() |
[17] |
Valiathan RR, Marco M, Leitinger B, et al. (2012) Discoidin domain receptor tyrosine kinases: new players in cancer progression. Cancer Metastasis Rev 31: 295-321. doi: 10.1007/s10555-012-9346-z
![]() |
[18] |
Ortega N, Werb Z (2002) New functional roles for non-collagenous domains of basement membrane collagens. J Cell Sci 115: 4201-4214. doi: 10.1242/jcs.00106
![]() |
[19] | Acerbi I, Cassereau L, Dean I, et al. Human breast cancer invasion and aggression correlates with ECM stiffening and immune cell infiltration. Integr Biol (Camb). [in press] |
[20] | Tao G, Levay AK, Peacock JD, et al. Collagen XIV is important for growth and structural integrity of the myocardium. J Mol Cell Cardiol 53: 626-638. |
[21] |
Mehner C, Radisky DC (2013) Triggering the landslide: The tumor-promotional effects of myofibroblasts. Exp Cell Res 319: 1657-1662. doi: 10.1016/j.yexcr.2013.03.015
![]() |
[22] |
Gehler S, Ponik SM, Riching KM, et al. (2013) Bi-directional signaling: extracellular matrix and integrin regulation of breast tumor progression. Crit Rev Eukaryot Gene Expr 23: 139-157. doi: 10.1615/CritRevEukarGeneExpr.2013006647
![]() |
[23] | Nistico P, Bissell MJ, Radisky DC (2012) Epithelial-mesenchymal transition: general principles and pathological relevance with special emphasis on the role of matrix metalloproteinases. Cold Spring Harb Perspect Biol 4. |
[24] |
Favreau AJ, Vary CP, Brooks PC, et al. (2014) Cryptic collagen IV promotes cell migration and adhesion in myeloid leukemia. Cancer Med 3: 265-272. doi: 10.1002/cam4.203
![]() |
[25] |
Emsley J, Knight CG, Farndale RW, et al. (2000) Structural basis of collagen recognition by integrin alpha2beta1. Cell 101: 47-56. doi: 10.1016/S0092-8674(00)80622-4
![]() |
[26] |
Espinosa Neira R, Salazar EP (2012) Native type IV collagen induces an epithelial to mesenchymal transition-like process in mammary epithelial cells MCF10A. Int J Biochem Cell Biol 44: 2194-2203. doi: 10.1016/j.biocel.2012.08.018
![]() |
[27] |
Leitinger B (2011) Transmembrane collagen receptors. Annu Rev Cell Dev Biol 27: 265-290. doi: 10.1146/annurev-cellbio-092910-154013
![]() |
[28] |
Noordeen NA, Carafoli F, Hohenester E, et al. (2006) A transmembrane leucine zipper is required for activation of the dimeric receptor tyrosine kinase DDR1. J Biol Chem 281: 22744-22751. doi: 10.1074/jbc.M603233200
![]() |
[29] |
Mihai C, Chotani M, Elton TS, et al. (2009) Mapping of DDR1 distribution and oligomerization on the cell surface by FRET microscopy. J Mol Biol 385: 432-445. doi: 10.1016/j.jmb.2008.10.067
![]() |
[30] |
Shrivastava A, Radziejewski C, Campbell E, et al. (1997) An orphan receptor tyrosine kinase family whose members serve as nonintegrin collagen receptors. Mol Cell 1: 25-34. doi: 10.1016/S1097-2765(00)80004-0
![]() |
[31] |
Vogel W, Gish GD, Alves F, et al. (1997) The discoidin domain receptor tyrosine kinases are activated by collagen. Mol Cell 1: 13-23. doi: 10.1016/S1097-2765(00)80003-9
![]() |
[32] |
Canning P, Tan L, Chu K, et al. (2014) Structural mechanisms determining inhibition of the collagen receptor DDR1 by selective and multi-targeted type II kinase inhibitors. J Mol Biol 426: 2457-2470. doi: 10.1016/j.jmb.2014.04.014
![]() |
[33] | Barker KT, Martindale JE, Mitchell PJ, et al. (1995) Expression patterns of the novel receptor-like tyrosine kinase, DDR, in human breast tumours. Oncogene 10: 569-575. |
[34] | Di Marco E, Cutuli N, Guerra L, et al. (1993) Molecular cloning of trkE, a novel trk-related putative tyrosine kinase receptor isolated from normal human keratinocytes and widely expressed by normal human tissues. J Biol Chem 268: 24290-24295. |
[35] | Karn T, Holtrich U, Brauninger A, et al. (1993) Structure, expression and chromosomal mapping of TKT from man and mouse: a new subclass of receptor tyrosine kinases with a factor VIII-like domain. Oncogene 8: 3433-3440. |
[36] |
Fu HL, Valiathan RR, Arkwright R, et al. (2013) Discoidin domain receptors: unique receptor tyrosine kinases in collagen-mediated signaling. J Biol Chem 288: 7430-7437. doi: 10.1074/jbc.R112.444158
![]() |
[37] |
Jin P, Zhang J, Sumariwalla PF, et al. (2008) Novel splice variants derived from the receptor tyrosine kinase superfamily are potential therapeutics for rheumatoid arthritis. Arthritis Res Ther 10: R73. doi: 10.1186/ar2447
![]() |
[38] |
Leitinger B (2014) Discoidin domain receptor functions in physiological and pathological conditions. Int Rev Cell Mol Biol 310: 39-87. doi: 10.1016/B978-0-12-800180-6.00002-5
![]() |
[39] |
Curat CA, Eck M, Dervillez X, et al. (2001) Mapping of epitopes in discoidin domain receptor 1 critical for collagen binding. J Biol Chem 276: 45952-45958. doi: 10.1074/jbc.M104360200
![]() |
[40] |
Abdulhussein R, McFadden C, Fuentes-Prior P, et al. (2004) Exploring the collagen-binding site of the DDR1 tyrosine kinase receptor. J Biol Chem 279: 31462-31470. doi: 10.1074/jbc.M400651200
![]() |
[41] | Leitinger B (2003) Molecular analysis of collagen binding by the human discoidin domain receptors, DDR1 and DDR2. Identification of collagen binding sites in DDR2. J Biol Chem 278: 16761-16769. |
[42] |
Ichikawa O, Osawa M, Nishida N, et al. (2007) Structural basis of the collagen-binding mode of discoidin domain receptor 2. EMBO J 26: 4168-4176. doi: 10.1038/sj.emboj.7601833
![]() |
[43] |
Xu H, Raynal N, Stathopoulos S, et al. (2011) Collagen binding specificity of the discoidin domain receptors: binding sites on collagens II and III and molecular determinants for collagen IV recognition by DDR1. Matrix Biol 30: 16-26. doi: 10.1016/j.matbio.2010.10.004
![]() |
[44] |
Carafoli F, Mayer MC, Shiraishi K, et al. (2012) Structure of the discoidin domain receptor 1 extracellular region bound to an inhibitory Fab fragment reveals features important for signaling. Structure 20: 688-697. doi: 10.1016/j.str.2012.02.011
![]() |
[45] |
Leitinger B, Kwan AP (2006) The discoidin domain receptor DDR2 is a receptor for type X collagen. Matrix Biol 25: 355-364. doi: 10.1016/j.matbio.2006.05.006
![]() |
[46] |
Konitsiotis AD, Raynal N, Bihan D, et al. (2008) Characterization of high affinity binding motifs for the discoidin domain receptor DDR2 in collagen. J Biol Chem 283: 6861-6868. doi: 10.1074/jbc.M709290200
![]() |
[47] |
Koo DH, McFadden C, Huang Y, et al. (2006) Pinpointing phosphotyrosine-dependent interactions downstream of the collagen receptor DDR1. FEBS Lett 580: 15-22. doi: 10.1016/j.febslet.2005.11.035
![]() |
[48] |
Wang CZ, Su HW, Hsu YC, et al. (2006) A discoidin domain receptor 1/SHP-2 signaling complex inhibits alpha2beta1-integrin-mediated signal transducers and activators of transcription 1/3 activation and cell migration. Mol Biol Cell 17: 2839-2852. doi: 10.1091/mbc.E05-11-1068
![]() |
[49] |
Lemeer S, Bluwstein A, Wu Z, et al. (2012) Phosphotyrosine mediated protein interactions of the discoidin domain receptor 1. J Proteomics 75: 3465-3477. doi: 10.1016/j.jprot.2011.10.007
![]() |
[50] |
Ruiz PA, Jarai G (2011) Collagen I induces discoidin domain receptor (DDR) 1 expression through DDR2 and a JAK2-ERK1/2-mediated mechanism in primary human lung fibroblasts. J Biol Chem 286: 12912-12923. doi: 10.1074/jbc.M110.143693
![]() |
[51] | L'Hote C G, Thomas PH, Ganesan TS (2002) Functional analysis of discoidin domain receptor 1: effect of adhesion on DDR1 phosphorylation. FASEB J 16: 234-236. |
[52] |
Huang Y, Arora P, McCulloch CA, et al. (2009) The collagen receptor DDR1 regulates cell spreading and motility by associating with myosin IIA. J Cell Sci 122: 1637-1646. doi: 10.1242/jcs.046219
![]() |
[53] |
Shintani Y, Fukumoto Y, Chaika N, et al. (2008) Collagen I-mediated up-regulation of N-cadherin requires cooperative signals from integrins and discoidin domain receptor 1. J Cell Biol 180: 1277-1289. doi: 10.1083/jcb.200708137
![]() |
[54] |
Hidalgo-Carcedo C, Hooper S, Chaudhry SI, et al. (2011) Collective cell migration requires suppression of actomyosin at cell-cell contacts mediated by DDR1 and the cell polarity regulators Par3 and Par6. Nat Cell Biol 13: 49-58. doi: 10.1038/ncb2133
![]() |
[55] |
Hansen C, Greengard P, Nairn AC, et al. (2006) Phosphorylation of DARPP-32 regulates breast cancer cell migration downstream of the receptor tyrosine kinase DDR1. Exp Cell Res 312: 4011-4018. doi: 10.1016/j.yexcr.2006.09.003
![]() |
[56] |
Hilton HN, Stanford PM, Harris J, et al. (2008) KIBRA interacts with discoidin domain receptor 1 to modulate collagen-induced signalling. Biochim Biophys Acta 1783: 383-393. doi: 10.1016/j.bbamcr.2007.12.007
![]() |
[57] | Dejmek J, Leandersson K, Manjer J, et al. (2005) Expression and signaling activity of Wnt-5a/discoidin domain receptor-1 and Syk plays distinct but decisive roles in breast cancer patient survival. Clin Cancer Res 11: 520-528. |
[58] |
Kim HG, Hwang SY, Aaronson SA, et al. (2011) DDR1 receptor tyrosine kinase promotes prosurvival pathway through Notch1 activation. J Biol Chem 286: 17672-17681. doi: 10.1074/jbc.M111.236612
![]() |
[59] |
Lu KK, Trcka D, Bendeck MP (2011) Collagen stimulates discoidin domain receptor 1-mediated migration of smooth muscle cells through Src. Cardiovasc Pathol 20: 71-76. doi: 10.1016/j.carpath.2009.12.006
![]() |
[60] |
Ongusaha PP, Kim JI, Fang L, et al. (2003) p53 induction and activation of DDR1 kinase counteract p53-mediated apoptosis and influence p53 regulation through a positive feedback loop. EMBO J 22: 1289-1301. doi: 10.1093/emboj/cdg129
![]() |
[61] |
Das S, Ongusaha PP, Yang YS, et al. (2006) Discoidin domain receptor 1 receptor tyrosine kinase induces cyclooxygenase-2 and promotes chemoresistance through nuclear factor-kappaB pathway activation. Cancer Res 66: 8123-8130. doi: 10.1158/0008-5472.CAN-06-1215
![]() |
[62] |
Ikeda K, Wang LH, Torres R, et al. (2002) Discoidin domain receptor 2 interacts with Src and Shc following its activation by type I collagen. J Biol Chem 277: 19206-19212. doi: 10.1074/jbc.M201078200
![]() |
[63] |
Yang K, Kim JH, Kim HJ, et al. (2005) Tyrosine 740 phosphorylation of discoidin domain receptor 2 by Src stimulates intramolecular autophosphorylation and Shc signaling complex formation. J Biol Chem 280: 39058-39066. doi: 10.1074/jbc.M506921200
![]() |
[64] |
Olaso E, Labrador JP, Wang L, et al. (2002) Discoidin domain receptor 2 regulates fibroblast proliferation and migration through the extracellular matrix in association with transcriptional activation of matrix metalloproteinase-2. J Biol Chem 277: 3606-3613. doi: 10.1074/jbc.M107571200
![]() |
[65] | Marcel V, Catez F, Diaz JJ (2015) p53, a translational regulator: contribution to its tumour-suppressor activity. Oncogene. [in press] |
[66] |
Petitjean A, Achatz MI, Borresen-Dale AL, et al. (2007) TP53 mutations in human cancers: functional selection and impact on cancer prognosis and outcomes. Oncogene 26: 2157-2165. doi: 10.1038/sj.onc.1210302
![]() |
[67] |
Vogel WF (2002) Ligand-induced shedding of discoidin domain receptor 1. FEBS Lett 514: 175-180. doi: 10.1016/S0014-5793(02)02360-8
![]() |
[68] |
Slack BE, Siniaia MS, Blusztajn JK (2006) Collagen type I selectively activates ectodomain shedding of the discoidin domain receptor 1: involvement of Src tyrosine kinase. J Cell Biochem 98: 672-684. doi: 10.1002/jcb.20812
![]() |
[69] |
Fu HL, Sohail A, Valiathan RR, et al. (2013) Shedding of discoidin domain receptor 1 by membrane-type matrix metalloproteinases. J Biol Chem 288: 12114-12129. doi: 10.1074/jbc.M112.409599
![]() |
[70] |
Shitomi Y, Thogersen IB, Ito N, et al. (2015) ADAM10 controls collagen signaling and cell migration on collagen by shedding the ectodomain of discoidin domain receptor 1 (DDR1). Mol Biol Cell 26: 659-673. doi: 10.1091/mbc.E14-10-1463
![]() |
[71] |
Huber MA, Kraut N, Beug H (2005) Molecular requirements for epithelial-mesenchymal transition during tumor progression. Curr Opin Cell Biol 17: 548-558. doi: 10.1016/j.ceb.2005.08.001
![]() |
[72] |
Knust E, Bossinger O (2002) Composition and formation of intercellular junctions in epithelial cells. Science 298: 1955-1959. doi: 10.1126/science.1072161
![]() |
[73] |
Yeh YC, Wu CC, Wang YK, et al. (2011) DDR1 triggers epithelial cell differentiation by promoting cell adhesion through stabilization of E-cadherin. Mol Biol Cell 22: 940-953. doi: 10.1091/mbc.E10-08-0678
![]() |
[74] |
Wang CZ, Yeh YC, Tang MJ (2009) DDR1/E-cadherin complex regulates the activation of DDR1 and cell spreading. Am J Physiol Cell Physiol 297: C419-429. doi: 10.1152/ajpcell.00101.2009
![]() |
[75] |
Yeh YC, Wang CZ, Tang MJ (2009) Discoidin domain receptor 1 activation suppresses alpha2beta1 integrin-dependent cell spreading through inhibition of Cdc42 activity. J Cell Physiol 218: 146-156. doi: 10.1002/jcp.21578
![]() |
[76] |
Xu H, Bihan D, Chang F, et al. (2012) Discoidin domain receptors promote alpha1beta1- and alpha2beta1-integrin mediated cell adhesion to collagen by enhancing integrin activation. PLoS One 7: e52209. doi: 10.1371/journal.pone.0052209
![]() |
[77] |
Imamichi Y, Menke A (2007) Signaling pathways involved in collagen-induced disruption of the E-cadherin complex during epithelial-mesenchymal transition. Cells Tissues Organs 185: 180-190. doi: 10.1159/000101319
![]() |
[78] |
Staudinger LA, Spano SJ, Lee W, et al. (2013) Interactions between the discoidin domain receptor 1 and beta1 integrin regulate attachment to collagen. Biol Open 2: 1148-1159. doi: 10.1242/bio.20135090
![]() |
[79] |
Rudra-Ganguly N, Lowe C, Mattie M, et al. (2014) Discoidin domain receptor 1 contributes to tumorigenesis through modulation of TGFBI expression. PLoS One 9: e111515. doi: 10.1371/journal.pone.0111515
![]() |
[80] |
Park CY, Min KN, Son JY, et al. (2014) An novel inhibitor of TGF-beta type I receptor, IN-1130, blocks breast cancer lung metastasis through inhibition of epithelial-mesenchymal transition. Cancer Lett 351: 72-80. doi: 10.1016/j.canlet.2014.05.006
![]() |
[81] |
Ozdamar B, Bose R, Barrios-Rodiles M, et al. (2005) Regulation of the polarity protein Par6 by TGFbeta receptors controls epithelial cell plasticity. Science 307: 1603-1609. doi: 10.1126/science.1105718
![]() |
[82] |
Buijs JT, Stayrook KR, Guise TA (2011) TGF-beta in the Bone Microenvironment: Role in Breast Cancer Metastases. Cancer Microenviron 4: 261-281. doi: 10.1007/s12307-011-0075-6
![]() |
[83] | de Jong JS, van Diest PJ, van der Valk P, et al. (1998) Expression of growth factors, growth-inhibiting factors, and their receptors in invasive breast cancer. II: Correlations with proliferation and angiogenesis. J Pathol 184: 53-57. |
[84] |
Walsh LA, Nawshad A, Medici D (2011) Discoidin domain receptor 2 is a critical regulator of epithelial-mesenchymal transition. Matrix Biol 30: 243-247. doi: 10.1016/j.matbio.2011.03.007
![]() |
[85] |
Xu J, Lu W, Zhang S, et al. (2014) Overexpression of DDR2 contributes to cell invasion and migration in head and neck squamous cell carcinoma. Cancer Biol Ther 15: 612-622. doi: 10.4161/cbt.28181
![]() |
[86] |
Maeyama M, Koga H, Selvendiran K, et al. (2008) Switching in discoid domain receptor expressions in SLUG-induced epithelial-mesenchymal transition. Cancer 113: 2823-2831. doi: 10.1002/cncr.23900
![]() |
[87] |
Siziopikou KP (2013) Ductal carcinoma in situ of the breast: current concepts and future directions. Arch Pathol Lab Med 137: 462-466. doi: 10.5858/arpa.2012-0078-RA
![]() |
[88] |
Turashvili G, Bouchal J, Ehrmann J, et al. (2007) Novel immunohistochemical markers for the differentiation of lobular and ductal invasive breast carcinomas. Biomed Pap Med Fac Univ Palacky Olomouc Czech Repub 151: 59-64. doi: 10.5507/bp.2007.010
![]() |
[89] |
Toy KA, Valiathan RR, Nunez F, et al. (2015) Tyrosine kinase discoidin domain receptors DDR1 and DDR2 are coordinately deregulated in triple-negative breast cancer. Breast Cancer Res Treat 150: 9-18. doi: 10.1007/s10549-015-3285-7
![]() |
[90] |
Vogel WF, Aszodi A, Alves F, et al. (2001) Discoidin domain receptor 1 tyrosine kinase has an essential role in mammary gland development. Mol Cell Biol 21: 2906-2917. doi: 10.1128/MCB.21.8.2906-2917.2001
![]() |
[91] |
Ren T, Zhang J, Liu X, et al. (2013) Increased expression of discoidin domain receptor 2 (DDR2): a novel independent prognostic marker of worse outcome in breast cancer patients. Med Oncol 30: 397. doi: 10.1007/s12032-012-0397-3
![]() |
[92] |
Morikawa A, Takeuchi T, Kito Y, et al. (2015) Expression of Beclin-1 in the Microenvironment of Invasive Ductal Carcinoma of the Breast: Correlation with Prognosis and the Cancer-Stromal Interaction. PLoS One 10: e0125762. doi: 10.1371/journal.pone.0125762
![]() |
[93] |
Zhang K, Corsa CA, Ponik SM, et al. (2013) The collagen receptor discoidin domain receptor 2 stabilizes SNAIL1 to facilitate breast cancer metastasis. Nat Cell Biol 15: 677-687. doi: 10.1038/ncb2743
![]() |
[94] | Lodillinsky C, Infante E, Guichard A, et al. (2015) p63/MT1-MMP axis is required for in situ to invasive transition in basal-like breast cancer. Oncogene. [in press] |
[95] |
Xiang S, Liu YM, Chen X, et al. (2015) ZEB1 Expression Is Correlated With Tumor Metastasis and Reduced Prognosis of Breast Carcinoma in Asian Patients. Cancer Invest 33: 225. doi: 10.3109/07357907.2015.1022258
![]() |
[96] |
Ling J, Kumar R (2012) Crosstalk between NFkB and glucocorticoid signaling: a potential target of breast cancer therapy. Cancer Lett 322: 119-126. doi: 10.1016/j.canlet.2012.02.033
![]() |
[97] |
Koh M, Woo Y, Valiathan RR, et al. (2015) Discoidin domain receptor 1 is a novel transcriptional target of ZEB1 in breast epithelial cells undergoing H-Ras-induced epithelial to mesenchymal transition. Int J Cancer 136: E508-520. doi: 10.1002/ijc.29154
![]() |
[98] |
Wolf K, Wu YI, Liu Y, et al. (2007) Multi-step pericellular proteolysis controls the transition from individual to collective cancer cell invasion. Nat Cell Biol 9: 893-904. doi: 10.1038/ncb1616
![]() |
[99] |
Neuhaus B, Buhren S, Bock B, et al. (2011) Migration inhibition of mammary epithelial cells by Syk is blocked in the presence of DDR1 receptors. Cell Mol Life Sci 68: 3757-3770. doi: 10.1007/s00018-011-0676-8
![]() |
[100] | Jonsson M, Andersson T (2001) Repression of Wnt-5a impairs DDR1 phosphorylation and modifies adhesion and migration of mammary cells. J Cell Sci 114: 2043-2053. |
[101] | Gavin BJ, McMahon AP (1992) Differential regulation of the Wnt gene family during pregnancy and lactation suggests a role in postnatal development of the mammary gland. Mol Cell Biol 12: 2418-2423. |
[102] |
Dejmek J, Dib K, Jonsson M, et al. (2003) Wnt-5a and G-protein signaling are required for collagen-induced DDR1 receptor activation and normal mammary cell adhesion. Int J Cancer 103: 344-351. doi: 10.1002/ijc.10752
![]() |
[103] |
Roarty K, Serra R (2007) Wnt5a is required for proper mammary gland development and TGF-beta-mediated inhibition of ductal growth. Development 134: 3929-3939. doi: 10.1242/dev.008250
![]() |
[104] |
Ribatti D, Crivellato E (2012) “Sprouting angiogenesis”, a reappraisal. Dev Biol 372: 157-165. doi: 10.1016/j.ydbio.2012.09.018
![]() |
[105] |
Artavanis-Tsakonas S, Rand MD, Lake RJ (1999) Notch signaling: cell fate control and signal integration in development. Science 284: 770-776. doi: 10.1126/science.284.5415.770
![]() |
[106] |
Dong Y, Li A, Wang J, et al. (2010) Synthetic lethality through combined Notch-epidermal growth factor receptor pathway inhibition in basal-like breast cancer. Cancer Res 70: 5465-5474. doi: 10.1158/0008-5472.CAN-10-0173
![]() |
[107] |
Thairu N, Kiriakidis S, Dawson P, et al. (2011) Angiogenesis as a therapeutic target in arthritis in 2011: learning the lessons of the colorectal cancer experience. Angiogenesis 14: 223-234. doi: 10.1007/s10456-011-9208-2
![]() |
[108] |
Yang Y, Sun M, Wang L, et al. (2013) HIFs, angiogenesis, and cancer. J Cell Biochem 114: 967-974. doi: 10.1002/jcb.24438
![]() |
[109] |
Ren T, Zhang W, Liu X, et al. (2014) Discoidin domain receptor 2 (DDR2) promotes breast cancer cell metastasis and the mechanism implicates epithelial-mesenchymal transition programme under hypoxia. J Pathol 234: 526-537. doi: 10.1002/path.4415
![]() |
[110] |
Rubinstein E (2011) The complexity of tetraspanins. Biochem Soc Trans 39: 501-505. doi: 10.1042/BST0390501
![]() |
[111] | Berditchevski F (2001) Complexes of tetraspanins with integrins: more than meets the eye. J Cell Sci 114: 4143-4151. |
[112] |
Powner D, Kopp PM, Monkley SJ, et al. (2011) Tetraspanin CD9 in cell migration. Biochem Soc Trans 39: 563-567. doi: 10.1042/BST0390563
![]() |
[113] | Miyake M, Nakano K, Ieki Y, et al. (1995) Motility related protein 1 (MRP-1/CD9) expression: inverse correlation with metastases in breast cancer. Cancer Res 55: 4127-4131. |
[114] |
Castro-Sanchez L, Soto-Guzman A, Navarro-Tito N, et al. (2010) Native type IV collagen induces cell migration through a CD9 and DDR1-dependent pathway in MDA-MB-231 breast cancer cells. Eur J Cell Biol 89: 843-852. doi: 10.1016/j.ejcb.2010.07.004
![]() |
[115] |
Hansen RK, Bissell MJ (2000) Tissue architecture and breast cancer: the role of extracellular matrix and steroid hormones. Endocr Relat Cancer 7: 95-113. doi: 10.1677/erc.0.0070095
![]() |