Analysis of spontaneous emergence of cell polarity with delayed negative feedback

1.   Introduction

Cell polarity refers to the asymmetric organization in the structure of cells and can enable cells to carry out specialized functions such as differentiation, migration and tissue development ^[1,2,3]. Disruption of cell polarity may lead to dysfunctionality of cells, which is usually a hallmark of human cancers ^[4,5]. Cell polarity always involves the localization of some specific signaling molecules to a proper location of the cell membrane ^[1,6], but the mechanism for achieving the localization remains controversial.

The budding yeast Saccharomyces cerevisiae has been developed as an attractive model system to study cell polarization as its generation time is short and useful experimental tools are available in this system ^[1]. The axis of polarized growth is regulated by budding location at where a newborn cell emerges from the original cell. The polarization event of budding yeast involves a key polarity protein, Cdc42 GTPase, which is highly conserved from yeasts to humans and is critical in polarity establishment ^[7,8]. Cdc42 localizes at a site of polarized growth on the plasma membrane and interacts with several proteins to trigger downstream processes, resulting in the emergence and growth of a bud ^[1]. In a wild-type cell, the Cdc42 localization generally occurs in response to spatial cues that is dependent on cell type. However, yeast cells can still bud in the absence of spatial cues, but the bud site is selected at a fully random manner, so-called symmetry breaking ^[9,10]. A lot of previous experimental works have studied polarization in the absence of Rsr1, which connects the spatial cue and the downstream polarization pathway ^[9,11,12,13].

Several mathematical models have been proposed to understand the mechanisms of achieving symmetry breaking in cell polarity ^{[9,14,15,16,17]}. More recently, in the geometric representation of a cell, studies ^[17,29] present bulk-surface reaction diffusion systems which incorporate cell membrane as surface and cytoplasm as bulk to investigate cell polarity in high dimensions. With the bulk-surface reaction diffusion models, classical linear stability analysis has been performed in literatures ^[17,30] and more recent approach, the local perturbation analysis, is also proposed ^[19,29]. Reaction-diffusion equation was always applied to model the budding process in yeast systems ^[14,15,16]. By the fact that the ratio of the diffusion rates of the cytoplasmic and membrane-bound Cdc42 is large, the Turing-type system became a possible mechanism for achieving symmetry breaking ^[14,16,18]. The key component of the Turing-type mechanism is the positive feedback loop in the cycle of Cdc42 and this concept was supported by a number of theoretical studies ^{[14,16,17,18,19]}.

Some experimental and theoretical studies indicated that negative feedback regulation exists during the formation of Cdc42 cluster ^{[20,21,22,23,24]}. Recent studies in ^[22,23,25] proposed that a delayed negative feedback may be important for maintaining the robustness of Cdc42 localization and the oscillating behavior of Cdc42 cluster during the process. However, cell polarization system with delayed negative feedback has not been studied with mathematical analysis in detail. In this paper, we first formulate a partial differential equation with positive and delayed negative feedback loops for studying cell polarization system. Then the Turing stability analysis is applied to identify the parameter conditions for achieving symmetry breaking during the process ^[14,26]. Our theoretical study provides us with several conditions for achieving a signaling cluster and simultaneously, these conditions include the constraint for the time delay in negative feedback loop. Numerical simulations are used to support our theoretical study and provide a full picture to understand the dynamics of cell polarization process. Also, our results support that the oscillating behavior of signaling cluster can be achieved by controlling the length of the time delay in negative feedback and the magnitude of positive feedback.

This paper is organized as follows. In Section 2, we present a reaction-diffusion model of cell polarization with positive and delayed negative feedbacks. In Section 3, we perform the Turing stability analysis to the model proposed in Section 2 to derive the conditions for which cell polarity may emerge. Section 4 contains the studies of our numerical simulations. Finally, the conclusion is presented in Section 5.

2.   Mathematical model for cell polarization with delayed negative feedback

Here we consider a partial differential equation describing the dynamics of a polarized signaling molecule on the cell membrane (Figure 1). A recent study ^[29] of bulk-surface reaction diffusion models suggested that the geometry of the cell may affect the dynamics of polarized signaling molecules but here we assume that the cell is a sphere which is consistent with the situation of budding yeast. Under the assumption applied in ^[9,22], we consider that all cytoplasmic signaling molecules are inactive and all membrane-bound signaling molecules are active. The variable we consider here is the membrane-bound active form of the signaling molecules and the cytoplasmic inactive form of this molecule is modeled implicitly through the conservation of total molecules. This type of model was well applied to study Cdc42-GTPase cycle in budding yeast ^[9,22]. Most of the GTPase cycles have a common feature that enables them to switch between their active and inactive forms. The activation process is usually initiated by hydrolysis and can be reversed by Guanine nucleotide exchange factors (GEFs), which cause the GDP to dissociate from the GTP. For the membrane-bound form, GDIs release the GDP from the cell membrane to the cytoplasm through binding to the GTPase and it can also be reversed through the action of GDI displacement factors.

The spatial domain in our model is the membrane of a cell of radius $R \ \mu m$ denoted by $M$, which is a sphere (two-dimensional domain), or, for simplicity, we can only consider the cross section of the cell membrane, which is a circle (one-dimensional domain) (shown in Figure 1A). We use a variable $a$ to represent the particle fractions of the membrane-bound active signaling molecules ^[9]. The dynamic of $a$ is governed by a reaction-diffusion equation with the feedback functions $F(\cdot, \cdot)$ and $G(\cdot)$:

with $\widehat{a}(t) = \int_M a(t, x) \ { dS_x}/ |M|$ and $\widehat{a^2}(t) = \int_M a^2(t, x) \ { dS_x}/ |M|$, respectively representing the average values of $a$ and $a^2$ over the cell membrane, and $|M|$ equals to the total area of the domain $M$. The first term of the right-hand side represents the diffusion of signaling molecules with the lateral surface diffusion rate coefficient $D_m$ and the Laplace operator $\nabla^2$ on the cell membrane.

In our model, we assume that the total number of signaling molecules in the whole cell is conserved. With the fact that $\hat{a}$ represents the total fractions of the membrane bound species, we obtain

where $\text{Frac}_c$ is the fraction of cytoplasmic inactive signaling molecules and $N$ is the total number of signaling molecules, including active and inactive forms, in the whole cell. Hence, by (2.2), $\text{Frac}_c = 1-\hat{a}$. Due to the fast cytoplasmic diffusion, we assume that signaling molecules are uniformly distributed throughout the cytoplasm. By assuming the activation rate is proportional to the fraction of cytoplasmic signaling molecules, the term $F(a, \widehat{a^2})(1-\hat{a})$ in (2.1) is applied to model the activation process in the system.

In budding yeast, there is a positive feedback loop to promote Cdc42 activation ^[9,22]. Here we assume that the activation rate is positively regulated by the active molecules $a$ and thus the function $F$ is an increasing function of $a$ (Figure 1B). In this paper, we define the feedback function as

This feedback function models multi-step cooperative interactions which has been used in several biological Turing type systems ^[26]. The nonlinear cooperativity is modeled by the term $a^2$, the degree of cooperativity is 2 ^[22]. The function form in (2.3) depends on a non-local term $\widehat{a^2}$ and the local density $a$. The feedback is mediated through a special type of molecules initially uniformly distributed in the cytoplasm, such as the Bem1 complex in the Cdc42 cycle of budding yeast ^[18,27]. The detailed derivation of this feedback function can be found in A.1.

The observed oscillation and fluctuation of Cdc42 cluster support that the delayed negative feedback is involved in the cell polarization system of budding yeast ^[22]. We assume that the deactivation function $G$ depends on the value of $a(t-\tau, x)$ with a delay time of $\tau$, as the deactivation rate varies with the activation level of Rga1, which may be regulated by the level of Cdc42 in budding yeast (Figure 1B). Here we apply a linear function to model the deactivation rate $G$ in (2.1):

The studies in ^[11] show that the necessary condition for achieving Cdc42 localization is that the rate of activation grows faster with the increase in active Cdc42 than the rate of deactivation, so that the linear function is used here instead of higher order functions.

3.   Linear stability analysis

In this section, the Turing stability analysis is applied to figure out the conditions of the parameters to achieve spontaneous cell polarization ^[26]. The analysis in this section can be applied for the system with general feedback functions $F(a, \widehat{a^2})$ and $G(a(t-\tau, x))$.

For studying the stability of a homogeneous steady state solution, we first study a homogeneous steady state solution $a_0$ of the system (2.1). Since $a_0$ is homogeneous over space, $\widehat{a_0^2} = a_0^2$ and the solution $a_0$ satisfies the following equation:

Since $F$ and $G$ are positive functions, the right-hand side of (3.1) is negative when $a_0 = 1$; the right-hand side of (3.1) is positive when $a_0 = 0$. By the intermediate value theorem, at least one homogeneous positive steady state solution $a_0$ exists between 0 and 1.

Here we define $a(t, x)$ as slightly perturbed functions from the homogeneous steady state $a_0$:

where the perturbation amplitude $\epsilon\ll1$ is much smaller than the value of $a_0$. After substituting (3.2) into Eq. (2.1) and applying Taylor expansion around $a_0$, the leading terms satisfy the following system:

where $a_1^* = a_1(t-\tau, x)$; $F_{X_1}$ and $F_{X_2}$ denote the partial derivatives of $F$ with respect to the first and the second arguments, respectively; $G_X$ denotes the first derivative of $G$. When the feedback functions (2.3) and (2.4) are considered, we obtain that $F_{X_1}$ is positive, $F_{X_2}$ is negative and $G_X$ is positive. A particular spatially periodic perturbation function, $a_1(t, x) = \alpha e^{\lambda t }E_w(x)$, is considered here. In the function, $\alpha$ is a nonzero parameters, $w$ is a non-negative integer and $E_w(x)$ is the $w$-th non-zero eigenfunction of the Laplace operator. Eq. (3.3) becomes

where

and the eigenvalue

where $R$ is the radius of the circle; $F_{X_1}, F_{X_2}, F, G, G_X$ are evaluated at $(a_0, a_0^2)$.

The emergence of cell polarity usually depends on the instability of the homogeneous steady state, which requires two conditions ^[16,26]:

(1) If the perturbation is spatially homogeneous ($w = 0$), the homogeneous steady state $a_0$ is linearly stable. This condition is equivalent to that the real parts of all eigenvalues $\lambda$ are negative when the wave number $w$ is zero. This condition ensures that a homogeneous solution starting from a constant initial condition close to $a_0$ will finally tend to $a_0$.

(2) The homogeneous steady state $a_0$ is linearly unstable under a perturbation with a positive wave number $w$. This condition is equivalent to that there exists a positive integer $w$ such that at least one $\lambda$ satisfying (3.4) has a positive real part.

These two conditions imply that the random perturbed homogeneous steady state are moving toward another inhomogeneous steady state for and only for positive wave lengths.

Eq. (3.4) has a form

where

and

According to ^[28], we can determine the signs of all possible $\lambda$ by the following theorem:

Theorem 1 (Theorem 4.7 from ^[28]). For Eq. (3.5), we have

(a) If $A_w+B > 0$, then there exists at least one $\lambda$ with positive real part;

(b) If $A_w+B < 0$ and $A_w-B\leq 0$, then all $\lambda$ have negative real parts;

(c) If $A_w+B < 0$ and $A_w-B > 0$, then there exists

such that (1) all $\lambda$ have negative real parts for $0 < \tau < \tau^*$, and (2) there exists at least one $\lambda$ with positive real part for $\tau > \tau^*$.

Condition (1)

By Theorem 1, the condition (1) is equivalent to that when $w = 0$, we have

or

where $A_0 = (F_{X_1} + 2a_0F_{X_2})(1-a_0)-G - F$ and $B = -G_Xa_0$.

Since $B$ is always negative, we have $A_0-B > A_0+B$ and $|B| = -B$ so

and

So the two situations for satisfying the condition (1) can be simplified as the following two cases:

or

From Cases (1a) and (1b), we can see that if the time delay $\tau$ is large and $A_0-B > 0$, the homogeneous steady state is not stable for a homogeneous perturbation. Under this situation, the steady state does not satisfy the condition (1) of Turing instability. For studying this situation, we consider the spatial homogeneous system from the model (2.1):

It is easy to show that the solution $\bar{a}$ is bounded between 0 and 1 with $\bar{a}(0) \in (0, 1)$. If all steady states of $\bar{a}$ are not stable for a large time delay $\tau$, the solution is under oscillation around some homogeneous steady states. Hence, for the original system (2.1), the instability for some positive wave number $w > 0$ may still contribute to achieve inhomogeneous pattern when the condition (2) is satisfied and the solution is close to the homogeneous steady state. Our simulation results shown in the later section will demonstrate that, with a large time delay, a cluster of signaling molecules can be formed from a homogeneous steady state with a small inhomogeneous perturbation even though all homogeneous steady state is not stable for $w = 0$.

Condition (2) By Theorem 1, the condition (2) is equivalent to the situation that for some positive integers $w > 0$, we have

or

where $A_w = -{ \sigma_w}D_m +F_{X_1}(1-a_0)-G$ and $B$ is defined above.

Since $A_w$ is decreasing with respect to $w$ for $w\geq 1$, the two cases can be simplified to the situation only for $w = 1$:

or

where $A_1 = - \sigma_1D_m+F_{X_1}(1-a_0)-G$.

To obtain Turing instability, at least one of Cases (2a) and (2b) has to be satisfied. The following two propositions provide the necessary conditions for obtaining at least one of Cases (2a) and (2b).

Proposition 2. For the system (2.1) with the feedback forms (2.3) and (2.4), Case (2a) is satisfied only if

Proof. Case (2a) can be written as

Since

Case (2a) implies $k_{on}k_2 > g_2$.

Proposition 3. For the system (2.1) with the feedback forms (2.3) and (2.4), Case (2a) or Case (2b) is satisfied only if

Proof. If Case (2a) or Case (2b) is true, we have $A_1-B > 0$ which can be written as

Since

the inequality (3.9) implies $k_{on}k_2 > \frac{1}{2}(g_1+ \sigma_1D_m)$.

The two propositions show that $k_{on}k_2$ has to be large enough for achieving Turing instability. Other than the necessary conditions, the following theorem provides a sufficient condition for determining a range of parameters in which Case (2a) is satisfied:

Theorem 4. Assume that $\sigma_1D_m < k_{on}k_1$. For the system (2.1) with the feedback forms (2.3) and (2.4), Case (2a) is satisfied if

and

The proof of Theorem 4 is based on the method used in ^[14] and the detailed proof is presented in A.2. It is worth noting that the necessary and sufficient conditions can be obtained by the formula for the roots of quadratic equation but the formula is too complicated and not applicable for parameter estimation.

The conditions obtained in Theorem 4 provide a good insight on the range of parameters for achieving Case (2a). For example, Theorem 4 implies that there exist a constant $k^* > 0$ such that Case (2a) is satisfied if $k_{on} > k^*$. Figure 2 displays the numerical results for determining how close between the conditions in Theorem 4 and the exact range for achieving Case (2a). In the figure, we consider 1000 sets of different $(k_2, g_2)$ generated uniformly within $[1,100]\times [1, 10]$, and choose the remaining parameters within the ranges: $\sigma_1D_m = [0.025, 0.225] /min$, $k_{on}\in[1, 5]/min$, $k_1 = 1$, $g_1 = 1/min$, which are based on some previous works ^[22,27].

Although the conditions in Theorem 4 (the brown regions in Figure 2A) are just sufficient conditions for Case (2a), they still have a good agreement with the exact range for Case (2a) (the brown regions in Figure 2B) and the difference appears only in a few sets of parameters (Figure 2C). The numerical results in Figure 2 show that the difference between the conditions in Theorem 4 and the exact region for Case (2a) can be minimized by reducing the membrane-bound diffusion coefficient $D_m$ which is usually small ($< 0.5 \mu m^2/min$) in budding yeast system ^[22,27].

4.   Numerical simulations

Linear stability analysis only focuses on the local behavior near the homogeneous steady state. Here we apply numerical simulations to study the long-term behavior of the system under different stability conditions studied in the previous section.

First, we use a computational simulation for one-dimensional model to study the ranges of the parameters for satisfying the two stability conditions (1) and (2). 10,000 sets of different $(k_2, g_2)$ are uniformly generated within the region $[1,100]\times [1,100]$. Other parameters are set as follows: $D_m = 0.1 \mu m^2/min$, $k_{on}\in[1, 5]/min$, $k_1 = 1$, $g_1 = 1/min$, $\tau\in[0.1, 1]min$. Here $D_m$ was obtained by the smallest value of $\sigma_1 D_m$ ($\sigma_1D_m = 0.025{/min}$ and $\sigma_1 = 1/R^2$) in Figure 2 which has given the conclusion that smaller $D_m$ can reduce the difference between conditions in Theorem 4 and the exact region for Case (2a). Figure 3 displays the stability diagrams for $k_{on} = 1{/min}, 3{/min}$ and $5{/min}$ and $\tau = 0.1{min}, 0.5{min}$ and $1{min}$. Figure 3A shows the ranges of the parameters that satisfy the condition (1) (the stability condition for homogeneous perturbation with $w = 0$); Figure 3B shows the ranges of the parameters that satisfy the condition (2) (the stability condition for inhomogeneous perturbation with $w = 1$).

According to the results in Figure 3, different combinations of parameters are chosen under different stability conditions. Here we fix the parameters $D_m = 0.1 \mu m^2/min$, $k_{on} = 3/min$, $k_1 = 1$ and $g_1 = 1/min$ ^[22,27]. For other parameters, we choose within the ranges $k_2\in[1,100]$, $g_2\in[1,100]{/min}$ and $\tau\in[0.1, 1]{min}$. Under this parameter setting, the system has a unique homogeneous steady state solution $a_0$. In all the simulations throughout this paper, the initial conditions are the homogeneous steady states with spatial perturbation ($\pm10\%$ perturbation):

where $\eta(x)$ is a function of uniformly distributed random number between -1 and 1. If the perturbation is large, the initial condition may be larger than 1. Within the ranges of the parameters we tested, we have verified that $1.1a_0 < 1$. In some cases, larger perturbation ($\pm20\%$) is also considered for further study and the results are consistent with what we observed here. We apply the second-order central difference approximation for the diffusion term, periodic boundary conditions for two end sides, Riemann sum for the definite integrals and the built-in function dde23 in the MATLAB for the temporal simulation with the constant time delay $\tau$. For the spatial discretization, the number of spatial points is 300 with uniform distribution and the radius of the circle is $R = 2\mu m$.

Figure 4 demonstrates the simulations for $g_2 = 5 {/min}$ which represents the cases that the magnitude of negative feedback is low. When $\tau = 0.1{min}$ (short time delay) and $k_2 = 20$ (weak positive feedback), the system does not satisfy the condition (2) and the numerical simulation in Figure 4A shows that the concentration of $a$ rapidly returns back to the homogeneous steady state. When the magnitude of positive feedback ($k_2$) increases to $50$ or higher, the system satisfies the two conditions (1) and (2) and inhomogeneous patterns can be always observed in the simulation (Figure 4A). Figure 4B shows that when $\tau$ is increased to be $0.5{min}$, a cluster of signaling molecules does not exist if the condition (2) is not satisfied, such as $k_2 = 20$ (weak positive feedback). When $k_2$ increases to $50$ or higher, a cluster of signaling molecules appears but is not stable at a certain location. In this case, the signaling cluster is traveling with constant speed and keeps moving until the end.

For stronger negative feedback ($g_2 = 50{/min}$ or $100{/min}$) with large enough time delay and positive feedback, a cluster of signaling molecules can appear but may not be stable at a certain location. Figure 5A demonstrates that when $\tau = 0.1{min}$ and the system does not satisfy the condition (2), the concentration of $a$ rapidly returns back to the homogeneous steady state even under a strong positive feedback.

Figures 5B-D show that how the dynamics of the signaling clusters change when the magnitude of the positive feedback increases from 5 to 50. Interestingly, when $k_2\geq15$, Case (2b) is satisfied and an inhomogeneous pattern is observed: for $k_2 = 15$, the concentration of a signaling cluster is oscillating during the process; when $k_2$ increases to 50, a traveling signaling cluster appears. In Figure 5D, when $k_2 = 50$, both Case (1a) and Case (1b) are not satisfied so the homogeneous steady state is unstable for a homogeneous perturbation $w = 0$. The simulation shows that a signaling cluster can be obtained in this case even though the homogeneous steady state is not stable when $w = 0$. But this situation is usually with larger time delay constant so the signaling cluster may not be stable at a certain location.

Figure 6 summarizes the simulation results for studying how the long-term behavior of the system depends on the values of $k_2$ and $\tau$ with two different levels of negative feedback, $g_2 = 5 {/min}, 50{/min}$. The results are consistent with our previous stability analysis. From this result, we observe that suitable ranges of the magnitude of positive feedback and time delay are important for achieving an oscillation of signaling cluster, observed in experiments ^[22]. It is worth noting that although the numerical results obtained here only focused on a simplified domain (the one-dimensional domain for the cross section of the cell), we can apply the results onto the the two-dimensional spherical surface. Figure 7 displays three examples of the simulations on the two-dimensional spherical surface. For the spatial discretization, the number of spatial points is 4098 with uniform distribution and the radius of the circle is $R = 2\mu m$. These three simulations are corresponding to the one-dimensional simulations shown in Fig 4A and Figures 5B and C. In Figure 7A, a stable cluster is formed until the end; in Figures 7B, C, the concentration of the signaling cluster is oscillating during the process. The results in these three cases are consistent with what observed in the one-dimensional simulations.

5.   Conclusion

Recent experimental studies demonstrated that delayed negative feedback regulation may play an important role in robustness of Cdc42 localization and the oscillating behavior of Cdc42 cluster during cell polarization ^[22]. However, detailed mathematical analysis is not well studied for this system.

In this paper, we have built a simple model of reaction-diffusion equation for cell polarization system which is regulated by positive and delayed negative feedback together. Our model involves general forms of positive and negative feedbacks. We have applied Turing stability analysis to analyze the conditions that can give rise to spontaneous cell polarization. Moreover, our numerical studies reveal that Cdc42 cluster can form but may not be spatially stable when the time delay $\tau$ is large. Also, our numerical results support that the oscillating behavior of Cdc42 cluster can be achieved by controlling the length of time delay and the magnitude of positive feedback.

The results in this paper provide parameter conditions for the existence of polarized solutions in the cell polarization system with delayed negative feedback. Furthermore, the analysis of the feedback can be easily extended to higher dimensional domain and provides insights to understand other similar biological systems in which cell polarity is established. Also, like the study in ^[22], our results can be used to study how the spatial cue level and the time delay affect the oscillating behavior by involving a spatial cue function $k_{cue}$ in the positive feedback term:

where

In the future work, our model can be combined with a moving boundary system for cell shape change during budding, then we can extend our study of Cdc42 localization to cell morphogenesis.

Acknowledgments

WCL was partially supported by a CityU Strategic Research Grant (Project No. 7004697) and a grant from the Research Grants Council of the Hong Kong Special Administrative Region, China (Project No. CityU 11303117).

Conflict of interest

All authors declare no conflicts of interest in this paper.

Supplementary

A.   Derivation of the positive feedback function and the proof of the theorem

A.1.   Derivation of the positive feedback function

We assume that the feedback molecules, as shown in Figure 1B, are initially uniformly distributed in the cytoplasm and then the active signaling molecules will recruit them to the cell membrane. On the other hand, the activation rate of the signaling molecules is proportional to the density of the membrane-bound feedback molecules. Therefore, we obtain the following equation for $c$:

where $(h_1+h_2 a^2)$ and $h_{off}$ are the recruitment rate and the dissociation rate of the feedback molecules, respectively; $(1-\hat{c})$ is the fraction of the cytoplasmic feedback molecules ($\hat{c}$ is the average value of $c$ over the membrane).

We know that the dynamics of the feedback molecules is much faster than that of the signaling molecules ^[18,27], the variable $c$ can be approximated by the solution of the quasi steady state equation of Eq. (A.1):

By integrating the above equation over the membrane, we can obtain the value of $\hat{c}$. After substituting $\hat{c}$ back into the equation, we have

where ${ k_{1}}$ and ${ k_{2}}$ equal to $h_1/h_{off}$ and $h_2/h_{off}$, respectively. We assume that the feedback strength is linearly proportional to the value of $c$, then we obtain the feedback function (2.3).

A.2.   Proof of the theorem

A.2.1.   Proofs of lemmas

Here we will state two lemmas, which will be used in Section A.2.2 for the proofs of Theorems 4. First, we define a function $f_{y}(a)$ used in the lemmas:

where $\gamma_1, \gamma_2, \gamma_4, \gamma_5, \gamma_6 > 0$, $\gamma_3\geq1$ and $0\leq D < \gamma_1\gamma_3$.

Lemma 5. Assume that $\gamma_2 > \gamma_6$ and $y\in[0, y_*]$ where $\gamma_2(1-\gamma_3y_* - \gamma_4y_*^2)-\gamma_6 = 0$. The function $f_{y}$ in (A.2) has the following properties:

1. $\min\limits_{a\geq 0}f_{y}(a)$ equals to

which is strictly decreasing with respect to $y$ for $y\in(0, y_*)$.

2. For each $y$, there exist at most two solutions in $\{a|a\geq 0\}$ satisfying $f_{y}(a) = 0$.

3. There exists a number $y_m$ in $[0, y_*)$ such that two smooth functions $a_1(y)$, $a_2(y)$ can be well defined in the domain $[y_m, y_*)$ and the following properties hold:

(a) $\min\limits_{a\geq 0}f_{y}(a)\leq 0$ for any $y\in[y_m, y_*)$;

(b) $f_y(a_1(y)) = f_y(a_2(y)) = 0$ for any $y\in[y_m, y_*)$;

(c) $a_1(y) > a_2(y)\geq 0$ for any $y\in(y_m, y_*)$;

(d) $a_1'(y) > 0$ and $a_2'(y) < 0$ for any $y\in(y_m, y_*)$;

(e) $\lim\limits_{y\rightarrow y_*}a_1(y) = \infty$ and $\lim\limits_{y\rightarrow y_*}a_2(y) = 0$;

(f) $\frac{df_y}{da}\left|_{a = a_1(y)}\right. > 0$ and $\frac{df_y}{da}\left|_{a = a_2(y)}\right. < 0$ for any $y\in(y_m, < /italic > < italic > y_*)$;

(g) if there is at least one solution in $a\geq 0$ for $f_0(a) = 0$, then $y_m = 0$;

(h) if there is no solution in $a\geq 0$ for $f_0(a) = 0$, then $a_1(y_m) = a_2(y_m)$, $\frac{df_{y_m}}{da}\left|_{a = a_1(y_m)}\right. = \frac{df_{y_m}}{da}\left|_{a = a_2(y_m)}\right. = 0$ and $\min\limits_{a\geq 0}f_{y_m}(a) = 0$.

Proof. 1. First we consider the first and second derivatives of $f_y$,

By Eq. (A.4), we show that $\frac{d^2f_{y}(a)}{da^2} > 0$ for $y\in[0, y_*)$ and the minimum point in $\{a|a\geq 0\}$ is at

with

By the given condition $D < \gamma_1\gamma_3$, we can show that $\min\limits_{a\geq 0}f_{y}(a)$ is strictly decreasing with respect to $y$.

2. Suppose that $y$ is a fixed number. If $\min\limits_{a\geq 0}f_{y}(a) > 0$, there is no solution $a\geq 0$ satisfying $f_y(a) = 0$.

If $\min\limits_{a\geq 0}f_{y}(a) = 0$, the minimum point

is one of the roots for $f_y(a)$. As $\frac{df_{y}(a)}{da} > 0$ for $a > \bar{a}$ and $\frac{df_{y}(a)}{da} < 0$ for $0\leq a < \bar{a}$, $f_y(a) > f_y(\bar{a})$ for any $a\neq \bar{a}$, and therefore $\bar{a}$ is the only solution of $f_y(a) = 0$.

If $\min\limits_{a\geq 0}f_{y}(a) < 0$, by the fact that $f_y(0) > 0$, $\lim\limits_{a\rightarrow \infty}f_y(a) > 0$ and the intermediate value theorem, we can show that there are at least two solutions satisfying $f_y(a) = 0$. As $\frac{df_{y}(a)}{da} > 0$ for $a > \bar{a}$ and $\frac{df_{y}(a)}{da} < 0$ for $0\leq a < \bar{a}$, $f_y(a) > f_y(\bar{a})$ for any $a\neq \bar{a}$. So there are only two roots of $f_y(a)$: one is in $\left[0, \bar{a}\right)$, and the other is in $\left(\bar{a}, \infty\right)$.

3. By the result of part 1, $\min\limits_{a\geq 0}f_{y}(a)$ tends to $-\infty$ as $y$ is close to $y_*$. If $\min\limits_{a\geq 0}f_{y}(a) > 0$ for $y = 0$, according to the intermediate value theorem, we can find $y_m$ such that $\min\limits_{a\geq 0}f_{y_m}(a)$ equals zero; if $\min\limits_{a\geq 0}f_{y}(a)\leq 0$ for $y = 0$, we define $y_m = 0$.

Since $\min\limits_{a\geq 0}f_{y}(a)$ is strictly decreasing with respect to $y$, and according to the results of part 2, $f_y(a) = 0$ has two solutions $a$ for any $y\in(y_m, y_*)$, so we can define two functions $a_1(y)$ and $a_2(y)$ that satisfy $f_y(a_1(y)) = f_y(a_2(y)) = 0$ and $a_1(y) > a_2(y)$ for any $y\in(y_m, y_*)$, that is,

The derivative of $f_y(a)$ with respect to $y$ is $-(\gamma_1+\gamma_2a^2)(\gamma_3+2\gamma_4y)+D$, which is always negative due to $\gamma_1\gamma_3 > D$, and $f_y(a)$ is a smooth function with respect to $y$ and $a$, and therefore we can apply the inverse function theorem to show that $a_1(y)$ and $a_2(y)$ are smooth functions. By the definitions and the proof of part 2, it is easy to verify the properties (a, b, c, f, g, h).

By property (b), we have $f_y(a_1(y)) = 0$ and $f_y(a_2(y)) = 0$. By differentiating these two equations with respect to $y$ on both sides, we have $-(\gamma_1+\gamma_2a_1(y)^2)(\gamma_3+2\gamma_4y)+D + \frac{df_y}{da}(a_1(y))a_1'(y) = 0$ and $-(\gamma_1+\gamma_2a_2(y)^2)(\gamma_3+2\gamma_4y)+D + \frac{df_y}{da}(a_2(y))a_2'(y) = 0$. Hence we obtain

By property (f) and $\gamma_1\gamma_3 > D$, we proved that $a_1'(y) > 0$ and $a_2'(y) < 0$, which completes the proof of property (d).

From the proof of part 2, we have $a_2\in\left[0, \bar{a}\right)$ and $a_1\in \left(\bar{a}, \infty\right)$. So we know that $a_1(y)$ tends to infinity as $y$ goes to $y_*$. Since $a = 0$ is the solution for $f_{y_*}(a) = 0$, we have $\lim\limits_{y\rightarrow y_*}a_2(y) = 0$, which completes the proof of property (e).

Lemma 6. Let $\Omega$ be a solution of

and assume that $\Omega > 0$ and $\gamma_2 > \gamma_6$. If

is satisfied, then $\frac{df_{a_0}}{da}\left|_{a = a_0}\right. > 0$ holds for any solution $a_0$ satisfying $f_{a_0}(a_0) = 0$.

For the proofs of Lemmas 6, we define two functions $S_1$, $S_2$ in the domain $[y_m, 1/\gamma_3)$:

where $a_1$, $a_2$ and $y_m$ are defined in Lemma 5.

Proof. There are two parts in the proof:

1. Prove that if $S_1(y_m) < 0$, $\frac{df_{a_0}}{da}\left|_{a = a_0}\right. > 0$ holds for any solution $a_0\geq 0$ satisfying $f_{a_0}(a_0) = 0$.

2. Prove that condition (A.6) implies $S_1(y_m) < 0$.

By combining these two results, we can prove that if the condition (A.6) is satisfied, then $\frac{df_{a_0}}{da}\left|_{a = a_0}\right. > 0$ holds for any solution $a_0\geq 0$ satisfying $f_{a_0}(a_0) = 0$.

Proof of part 1: Suppose that $S_1(y_m) < 0$. Since $a_1(y)\geq a_2(y)$, we get $S_2(y_m)\leq S_1(y_m) < 0$. By $a_2'(y) < 0$ (Lemma 5(3c)), we have $S_2' < 0$, which means that $S_2$ is a decreasing function. Since $S_2(y_m) < 0$ and $S_2$ is a decreasing function, $S_2(y) < 0$ for all $y\in[y_m, y_*)$, and there is no solution to $S_2(y) = 0$.

According to Lemma 5 and the definitions of $S_1$ and $S_2$, all solutions $a_0\geq 0$ for $f_{a_0}(a_0) = 0$ have to satisfy $S_1(a_0) = 0$ or $S_2(a_0) = 0$. Since $S_1(y_m) < 0$ implies that there is no solution satisfying $S_2(y) = 0$, all solutions $a_0\geq0$ for $f_{a_0}(a_0)$ have to satisfy $S_1(a_0) = 0$ and therefore $\frac{df_{a_0}}{da}\left|_{a = a_0}\right. > 0$ according to Lemma 5(3f).

Proof of part 2: Suppose that condition (A.6) is satisfied, by Lemma 5(1), we have

If $0 < y < \Omega$, we have

Lemma 5(3a) implies that $\min\limits_{a\geq 0}f_{y_m}(a) > 0$ so $y_m$ is larger than $\Omega > 0$. Then we apply Lemma 5(3h) to show that there is no solution with $a\geq 0$ such that $f_0(a) = 0$.

By Lemma 5(3b, h), we know that ($y_m, a_1(y_m))$ satisfies the following two equations:

where $a_m = a_1(y_m)$.

After multiplying (A.8) and (A.9) by $n$ and $a_m$, respectively, we have

Subtracting (A.10) by (A.11), we obtain

which leads to

By substituting (A.12) into $S_1(y_m)$, we obtain

By applying $y_m > \Omega > 0$, $D < \gamma_1\gamma_3$ and condition (A.6),

which, coupled with (A.13), implies that $S_1(y_m) < 0$.

A.2.2.   Proof of Theorem 4

Proof. First, we set $\gamma_1 = k_{on}k_1$, $\gamma_2 = k_{on}k_2$, $\gamma_3 = \frac{k_{on}+g_1}{k_{on}} > 1$, $\gamma_4 = \frac{g_2}{k_{on}}$, $\gamma_5 = g_1$, $\gamma_6 = g_2$ and $D = \sigma_1D_m$ used in the lemmas.

Solving (A.5) in Lemma 6, we obtain

By the condition (3.10) of Theorem 4

we have $\Omega > 0$ and $\gamma_2 = k_{on}k_2 > \gamma_6 = g_2$. Also the other condition (3.11) of Theorem 4

implies the condition (A.6). By Lemma 6, we have

holds for any $a_0$ satisfying

Eq. (A.15) implies that

Now we get that $a_0$ is a homogeneous steady state solution for $a$ in system (2.1) with the feedback forms (2.3) and (2.4) if and only if $a_0$ satisfies (A.15).

Also, (A.16) and (A.15) imply that

and then combining (A.14), we have

which is equivalent to Case (2a).