## Abstract

Cell polarization is a ubiquitous process which results in cellular constituents being organized into discrete intracellular spatial domains. It occurs in a variety of cell types, including epithelial cells, immune system cells and neurons. A key player in this process is the Par protein family whose asymmetric localization to anterior and posterior parts of the cell is crucial for proper division and cell fate specification. In this paper, we explore a stochastic analogue of the temporal model of Par protein interactions first developed in Dawes & Munro (Dawes and Munro 2011 *Biophys. J.* 101, 1412–1422. (doi:10.1016/j.bpj.2011.07.030)). We focus on how protein abundance influences the behaviour of both the deterministic and stochastic versions of the model. In Dawes & Munro (2011), it was found that bistable behaviour in the temporal model of Par protein led to the existence of complementary domains in the corresponding spatio-temporal model. Here, we find that the corresponding temporal stochastic model permits switching behaviour (the model solution ‘jumps’ between steady states) for lower protein abundances, whereas for higher protein abundances the stochastic and deterministic models are in good agreement (the model solution evolves to one of two steady states). This led us to the testable hypothesis that cells with lower abundances of Par protein may be more sensitive to external cues, whereas cells with higher abundances of Par protein may be less sensitive to external cues. In order to gain more control over the precise abundance of Par protein, we proposed and explored a second model (again, examining both deterministic and stochastic versions) in which the total number of Par molecules is conserved. We found that this model required an additional dimerization reaction in the cytoplasm in order for bistable and switching behaviour to be found. Once this additional reaction was included, we found that both the first and second models gave qualitatively similar results but in different regions of the parameter space, suggesting a further regulatory mechanism that cells could potentially use to modulate their response to external signals.

## 1. Introduction

Intracellular polarization, where cells asymmetrically localize specific factors to opposite ends of the cell in order to specify a spatial axis, is a common and fundamental biological phenomenon. The spatial organization of these domains is crucial for a wide range of processes including directed motility, embryonic development and epithelial tissue integrity [1,2]. The establishment of these biochemically and dynamically distinct domains relies on core modules of polarity determinant proteins which are highly conserved across eukaryotic organisms [2,3].

In this study, we focus on the Par proteins, a key family of polarity determinants first characterized in the nematode worm *Caenorhabditis elegans* [4,5]. In early embryos of *C. elegans*, polarization occurs shortly after fertilization, with PAR-3, PAR-6 and the atypical protein kinase (aPKC) localizing to the anterior half of the embryo, and PAR-1, PAR-2 and tumour-suppressor protein LGL localizing to the posterior half, shown in figure 1. The protein PAR-3 contains a homo-oligomerization domain, allowing the anterior Par proteins to form higher order complexes [7,8]. In their active unphosphorylated form, the Par proteins bind to the cortex, an actin-based structure located on the periphery of the cell. The anterior and posterior Par proteins interact antagonistically via mutual phosphorylation, causing the phosphorylated protein to fall off the cortex. The inactive phosphorylated Par proteins then freely diffuse in the cytoplasm. The asymmetric domains in the early *C. elegans* embryo, which are maintained for approximately 10 min as the cell prepares for first division, are critical for continued embryonic development [9]. A similar asymmetric localization is observed in other polarizing cell types [1,2]. As these protein asymmetries are crucial for the proper functioning of polarized cells, a key challenge is delineating how the organization of Par protein asymmetries responds to external signals.

Intracellular polarization is observed in a wide variety of biological contexts, and the range of biological processes that rely on polarization is equally varied. In some contexts, polarized cells remain highly sensitive to additional signals that can disrupt the spatial localization, while other polarized cell types are less sensitive to disruptive signals. We refer to these two modes of polarization as ‘soft’ and ‘hard’.

Soft polarization refers to cells that remain sensitive to external signals and can reorient their polarized domains in response to that signal. This is commonly seen in motile cells such as *Dictyostelium* that will lose polarization, transiently becoming unpolarized, then repolarize in the direction specified by the new signal, which can be either chemical or mechanical [10,11]. Similarly, neutrophils, a component of the innate immune system, will reorient their polarized domains in response to a new signal [12,13].

By contrast, hard polarization refers to those polarized cells that maintain their polarized domains even in the face of chemical or mechanical perturbations. For instance, early *C. elegans* embryos are able to maintain stable polarized domains even after complete loss of the actin cytoskeleton [14] or loss of microtubules nucleated from centrosomes [15]. Epithelial cells also exhibit hard polarization; loss of epithelial polarity has grave consequences for tissue integrity and is a hallmark of metastatic cancer [16].

A large number of deterministic mathematical models been proposed to study polarization, particularly in the context of cell motility. These models, ranging from biochemically based to Turing-type models, are reviewed and compared in an excellent recent article (see [17]). These models have made a number of predictions about key interactions and molecules needed for polarization. Analysis of stochastic models has provided insights into the minimum number of molecules needed to initiate polarization [18,19].

A smaller number of mathematical models have been developed to explicitly study Par protein dynamics in polarization. These Par protein models explored initiation and maintenance of polarization in the early *C. elegans* embryo using reaction–diffusion models [14,20]. Previous models proposed by our group have focused on the establishment and maintenance of Par protein domains in the early embryo [6,21]. These previous Par protein studies have shed light on critical interactions between Par proteins and mechanical structures in the cell, notably the actin cytoskeleton. While some previous investigations, particularly when considering cell motility, have demonstrated whether a model cell is able to change polarity orientation in response to additional signals [17,22], there has not yet been an in-depth exploration of the conditions under which polarity orientation can be switched. In the context of this investigation, we use ‘switching’ to mean transient depolarization followed by re-establishment of polarized domains with a new axis of orientation.

In this paper, we build on the model of Par protein dynamics proposed in [6] to study polarization dynamics. That model consisted of a system of ordinary differential equations (ODEs) that assumed cytoplasmic protein levels were at quasi-steady state (QSS model). In that paper, a strictly deterministic version of the model was analysed to determine the conditions under which polarized domains could be stably maintained. Here, we build on that previous work by analysing the QSS model in both deterministic and stochastic implementations to investigate the conditions under which the orientation of polarized domains can be switched, particularly in response to persistent noise in the system. We further propose an extension of this model that assumes the total amount of protein (cytoplasmic plus cortical) is conserved (conservation model). We analyse and simulate both versions of the Par protein model, the QSS and conservation models, using deterministic and stochastic implementations. When we vary the size of the cytoplasmic pool in the QSS model or the amount of total protein in the conservation model, we find the deterministic formulations of both models are able to support bistable solutions for a range of values. These bistable solutions correspond to complementary protein distributions observed in the early embryo and are consistent with our earlier findings [6]. However, stochastic implementations of the QSS model and the conservation model are able to produce switching behaviours for intermediate protein levels, meaning that the system can alternate between two different protein configurations as a result of noise in the system. This suggests that a cell may be able to tune its sensitivity to external signals by modulating the amount of protein in the cell, and we propose experimentally testable hypotheses to investigate this behaviour *in vivo*.

## 2. Constant cytoplasmic pool model

In this section, we introduce our modelling assumptions and the Par protein interactions we consider. We also present the corresponding deterministic ODE model originally presented in [6] which we refer to as ‘R1’ for the rest of the paper. The parameters used are listed in table 1.

### 2.1. Model reactions

We keep track of three different forms of anterior Par proteins and one form of posterior Par protein. The model variables are as follows: *A*_{1}, cortical monomeric ParA; *A*_{10}, cortical dimeric ParA (singly bound); *A*_{11}, cortical dimeric ParA (doubly bound) and P, cortical ParP. In other words, *A*_{10} is a ParA dimer with only one molecule bound to the cortex, while *A*_{11} is a ParA dimer with both molecules bound to the cortex. We assume mass action kinetics for all reactions in the model. We assume the following quantities are at QSS and so are included in the model but do not vary over time: *A _{y}*, cytoplasmic ParA monomers;

*A*

_{2y}, cytoplasmic ParA dimers;

*P*, cytoplasmic ParP. All the interactions are summarized schematically in figure 2.

_{y}To begin, we assume that cytoplasmic monomeric ParA can bind to and dissociate from the cortex 2.1

Monomeric cortical ParA can then undergo a homodimerization reaction, creating a dimer with both molecules bound to the cortex 2.2

We also assume that monomeric cytoplasmic ParA can form a dimer with cortical ParA, creating a dimer with only one molecule bound to the cortex 2.3

The monomeric cortical ParA and both forms of dimeric ParA can be phosphorylated by ParP, at a rate *r _{A}* (after which we no longer keep track of their dynamical evolution)
2.4
2.5
2.6We also assume that cytoplasmic ParA dimers can bind to and dissociate from the cortex, in a sequential manner, i.e. one molecule at a time, which is encapsulated by the following reactions:
2.7and
2.8In a similar manner to anterior Par proteins, we assume that cytoplasmic monomeric ParP can bind to and dissociate from the cortex
2.9

We also assume that posterior proteins can be phosphorylated (and thus removed from the system) by all forms of the anterior protein module
2.10
2.11
2.12We consider the same parameter set that was explored in R1 (see column 2 of table 1). We also consider two additional parameter sets in which we increase *A _{y}* and

*P*(table 1). By increasing the abundance of proteins contained in the cytoplasmic pools,

_{y}*A*and

_{y}*P*, we effectively study the influence of copy number (protein abundance) on the behaviour of the model. We consider three parameter sets which we loosely classify as corresponding to ‘low-’, ‘medium-’ and ‘high-’ protein copy numbers. To correlate protein number to concentration, we estimate the volume of the

_{y}*C. elegans*embryo using the formula for an ellipsoid. It is known that the length of the embryo is approximately 60 µm, while the diameter is approximately 30 µm, making the average cell volume of a

*C. elegans*embryo approximately 60 µm or 6 × 10

^{−14}l [23,24].

### 2.2. Mean-field ordinary differential equation approach

Using the law of mass action, the corresponding deterministic ODE system can be derived from the reactions we presented in §2.1 and written as follows:
2.13
2.14
2.15
2.16which are subject to the following initial conditions:
where *A*_{0}, *A*_{10}, , and *P*_{0} are positive constants. In order to compare concentrations (real numbers) with discrete numbers of molecules (integers), we need to define how to convert concentrations to numbers. Concentrations are usually measured in M (moles l^{−1}). Avogadro's constant, , gives the number of molecules in a mole. Hence a concentration, for example *a*_{1}(*t*)M (the monomer form of anterior Par protein), in a fixed volume of vol litres corresponds to *a*_{1}(*t*)*n*_{A}vol molecules of that species. For the remainder of the paper, we present results in terms of numbers of molecules. In all figures, we plot total ParA and total ParP numbers (where total ParA is the sum of the three different ParA forms).

### 2.3. Linear stability analysis

Steady states of our ODE system can be written by setting the left-hand side of equations (2.13)–(2.16) equal to 0 and solving the resulting algebraic system. To accomplish this, we make use of the numeric algebraic solver ‘vpasolve’ in Matlab. For the parameter sets presented in table 1, it is straightforward to show there are three biologically feasible steady states for each parameter set and seven steady states in total. For example, using parameter set 1, the biologically feasible steady states are

The first steady state corresponds to the case where there is high ParP and low total ParA, which is similar to the protein distribution observed at the posterior pole of a polarized embryo. The second steady state corresponds to the case where there is low ParP and high total ParA, which is similar to the protein distribution observed at the anterior pole of a polarized embryo. The last steady state represents a state where there are relatively low abundances of all four protein species. By linearizing about each steady state, we obtain the following linear system:

where *A _{i}* is the Jacobian matrix of the system evaluated at the corresponding steady state. The eigenvalues of

*A*determine the stability of each steady state. It can be shown that the first two steady states are stable and the third is unstable. The system can also be described as possessing two distinct basins of attraction, being separated by a separatrix which passes through a saddle point (the third steady state). We find similar behaviour for the two other parameter sets in table 1, i.e. the linear stability analysis yields two stable steady states (with one steady state producing high Par P, low total ParA and vice versa for the other steady state) and a third unstable steady state.

_{i}In figure 3, we show the result of exploring the (*A _{y}*,

*P*) parameter space on the steady-state solutions and their stability. For the majority of the parameter space considered, there are three biologically feasible steady states (with two stable and one unstable). However, for certain regions, there exists only a single biologically feasible steady state (e.g. high

_{y}*P*and low

_{y}*A*) and for other parts there are no biologically feasible steady states (e.g. high

_{y}*A*and low

_{y}*P*). The white space in the subfigures indicates regions where no biologically feasible solution exists (i.e. the steady state is negative or complex). With respect to the first steady state, as

_{y}*A*is increased the abundance of total ParA increases (provided

_{y}*P*is below the value 3) and the abundance of ParP decreases. Whereas if

_{y}*P*is increased the abundance of ParP increases and total ParA decreases. For the second steady state, similar behaviour is observed but the third steady state exhibits counterintuitive behaviour. For the third steady state, as

_{y}*P*is increased total ParA levels increase, and as

_{y}*A*is increased ParP levels increase.

_{y}Figure 4 shows how the number of biologically feasible steady states changes with respect to variation in the sizes of the cytoplasmic pools, *A _{y}* and

*P*. For most of the parameter space shown, there are three steady states but we note that for large

_{y}*A*values relative to

_{y}*P*and vice versa, only one steady state can exist. This suggests that the cytoplasmic pools have to be roughly the same size in order to observe coexistence of complementary domains.

_{y}### 2.4. Numerical simulations of deterministic and stochastic models

In this section, we present solutions of the ODE model and compare them with solutions from the stochastic model. Equations (2.13)–(2.16) were solved numerically using the Matlab stiff ODE solver ode15 s and realizations of the stochastic model were computed using Gillespie's algorithm [25]. We stress that it is only the protein abundance that we alter in simulations, i.e. values of *A _{y}* and

*P*. The initial conditions we choose take the following form: 2.17 2.18 2.19 2.20These initial conditions are chosen arbitrarily to show how the solutions evolve for different initial quantities of the model species. The same initial conditions are used for both deterministic and stochastic models.

_{y}#### 2.4.1. Stochastic model does not produce bistability for low-protein levels

For parameter set 1, we find bistable behaviour in the ODE model which is consistent with the findings of R1 and the steady-state exploration we presented in §2.3. That is, depending on the initial condition, the solution tends to one of two steady states. This is shown clearly in the first column of the first row of figure 5, where we present a phase plane showing the result of long-time numerical simulations stemming from four different initial conditions. The ODE is solved for a time corresponding to 10 000 s (all other simulations are solved for this same time). The black dots represent the solution value at the final time step—in other words, they represent the state the solution converges to. Two possible steady states exist, one where ParP is high and total ParA is low and another where total ParA is high and ParP is low. However, the solution yielded by the corresponding stochastic model is not in agreement. The plot presented in the second column of the first row of figure 5 shows the evolution of the cumulative mean (or moving average) of total ParA and ParP numbers computed from a long-time simulation using the Gillespie algorithm. We define the cumulative mean as follows:
where *t* is the current time in the stochastic simulation, *τ* represents the time step and *x* represents the current value of a model variable. The cumulative mean is useful for uncovering general trends in the time series, smoothing out short-term fluctuations and finding out which mean value a time series converges to. It is clear that the stochastic model can only attain a single steady state which appears to have approximately equal amounts of total ParA and ParP (on average) regardless of initial conditions. We performed 1000 additional long-time simulations for random initial conditions to confirm this. In a low molecule number regime, there is a qualitative difference between the solutions of the stochastic and deterministic cases, potentially due to the high amount of noise from the low-protein abundance which traps the system in this state.

#### 2.4.2. Stochastic model produces switching behaviour for intermediate protein levels

For parameter set 2 (see column 3 of table 1), the long-time numerical solution of the ODE model again agrees with steady-state values shown in figure 3. There is clearly bistable behaviour exhibited in the ODE model and the particular steady state obtained once again depends on the initial condition chosen (see first column, second row of figure 5). Interestingly, we find very different behaviour in the corresponding stochastic model. In the second column, second row of figure 5, we can see that the cumulative mean for the stochastic simulation tends to the high total Par A, low ParP steady state for some initial conditions while for other initial conditions it tends to a value somewhere between the two different steady states of the corresponding ODE model. By inspecting the time series of the third initial condition (shown in green) which converges to a value between the two steady states, we find that the solutions ‘jumps’ between the two different steady states of the corresponding ODE model. This switching behaviour is due to the level of intrinsic noise being enough to continually perturb the stochastic solution from staying in one of the two steady states. We note that we repeated the simulations that produced the second column, second row of figure 5 several times and each time found the cumulative mean to converge to a different mean which was located either at one of the two steady states or some value along the straight line connecting the two steady states. Hence, for ‘medium’ levels of proteins, there is a qualitative difference in the solutions of the stochastic and deterministic models. Finally, in the third column, second row of figure 5, we present one clear example of a trajectory which undergoes a stochastic ‘jump’ at approximately *t* = 600 s causing the steady state to change from high ParP, low ParA to high ParA, low ParP.

#### 2.4.3. Stochastic and deterministic models are in good agreement for large abundances of proteins

Finally, for parameter set 3 (column 4 of table 1) which corresponds to an embryo with a relatively high abundance of ParA and ParP molecules, we observe similar behaviour for both the stochastic and deterministic cases (see third row of figure 5). The solutions of the stochastic model now appear to depend strongly on the initial conditions which are in good agreement with the corresponding deterministic model. This is consistent with other studies that have shown stochastic systems agree well with their deterministic counterparts as the number of molecules increases [26,27]. No other states are occupied except those that are near the initial conditions and noise is unable to ‘kick’ the solution into another steady state. That is, we do not observe the same switching behaviour that we observed for moderate levels of protein.

## 3. Conservation of total protein model

In this section, we deviate from the model of R1 to consider a model in which the cytoplasmic pools are not at QSS, but rather we assume the total amount of protein is conserved. In other words, we no longer assume that the cytoplasmic pools are constant, but instead we assume they change with respect to time and the total amount of protein at any moment in time is constant and is equal to the sum of the initial quantities of the model species where is the initial number of anterior Par monomers, is the initial number of anterior Par dimers and is the initial number of posterior Par dimers. By reformulating the model in this way, we have greater control over the total abundance of protein and we can attribute model behaviours to protein abundances more precisely. For the case of the ODE model, we append three additional ODEs to the system given by equations (2.13)–(2.16) to yield the following system: 3.1 3.2 3.3 3.4 3.5 3.6 3.7which are subject to the following initial conditions:

In the Gillespie algorithm, in order to account for conservation of total protein numbers, we update the reaction propensities accordingly. For example, if reaction (2.1) fired, we would decrease the cytoplasmic pool, *A _{y}*, by one (making the reaction less likely to occur in the next time step), whereas for the stochastic model considered in §2, it would be left unchanged. To make this clear, we rewrite reactions (2.1), (2.4)–(2.7) and (2.9)–(2.12) as follows:
3.8
3.9
3.10
3.11
3.12
3.13
3.14
3.15
3.16
3.17

### 3.1. Linear stability analysis of protein conservation model

As we did in §2.2.1, we present here the steady states of the protein conservation model by setting the left-hand side of equations (3.1)–(3.7) equal to 0 and solving the resulting algebraic system (again using Matlab's ‘vpasolve’ function). By making use of the conservation of protein, we can define the following three variables (the total amounts of ParA monomers, ParA dimers and ParP proteins) which we use to simplify the steady-state system which allows us to rewrite the cytoplasmic pools as

Therefore, we can write the steady-state system as
3.18
3.19
3.20
3.21
3.22Surprisingly, the conservation of protein model appears to prevent the appearance of bistable solutions in the deterministic model (figure 6). There is only one biologically feasible steady state which is unstable for large values of *P _{T}*. This suggests that the existence of constant cytoplasmic pools is needed in order for bistable solutions to exist in the deterministic model. The solution of the conservation of the protein model always tends to a low ParA, high ParP state. We confirm this behaviour by looking at the numerical solution of the ODE model in §3.2 as well as exploring the corresponding stochastic model.

### 3.2. Numerical simulations of deterministic and stochastic protein conservation models

We use the same methods of numerical simulation as we did for the constant cytoplasmic pool model. However, this time we choose initial conditions that take the following form:
3.23
3.24
3.25
3.26where *A*_{0} and *P*_{0} are the initial numbers of ParA proteins and ParP proteins, respectively. Note that we set the ratio of posterior to anterior Par Protein abundances to be 3 : 10, in other words, for each initial condition there are three posterior Par proteins for every 10 anterior Par Proteins. This was motivated by figure 9, which revealed that three steady states are only possible when the abundance of ParP is far less than ParA. As for the constant cytoplasmic pool model, we use the same initial conditions for the both deterministic and stochastic protein conservation models. Again, we explore three different parameter regimes which correspond to ‘low’, ‘medium’ and ‘high’ total numbers of protein in an embryo. Specifically, we choose *A*_{0} and *P*_{0} to be equal to 1, 15 and 150, which correspond to total protein numbers equal to 13, 195 and 1950. Figure 7 shows the results of the numerical simulations for the deterministic model (first column) and stochastic model (second and third columns). In contrast to the constant cytoplasmic pool model, we find the two different modelling approaches are in good agreement regardless of the abundance of protein. As we found in our steady-state analysis of the deterministic system, there appears to be one steady state regardless of the parameter set—namely, high ParP, low ParA. We note that in the low-protein copy number regime (first row) that the total ParA number is not zero and in fact is in good agreement with the constant cytoplasmic pool model (compare the first row of figures 7 and 5).

### 3.3. Linear stability analysis of protein conservation model with cytoplasmic dimerization

Here, we adapt the protein conservation model to introduce two additional reactions in the cytoplasm: the dimerization of ParA monomers and the separation of ParA dimers into monomers. We achieve this by adding two additional reactions to equations (3.5) and (3.6), which are changed to become
3.27and
3.28We also account for dimerization in the cytoplasmic pools for the stochastic model, which is captured by the following reaction:
Note that we assume for the sake of simplicity that the dimerization rate in the cytoplasm is the same as the dimerization rate for cortex-bound species. The addition of dimerization to the ParA cytoplasmic pools dramatically changes the solution landscape. In figure 8, we present the result of exploring the (*A*_{T}, *P*_{T}) parameter space on the steady-state solutions and their stability. As in the case of the original model with cytoplasmic pools (see §2.3), we find that there can be either one or three biologically feasible steady states depending on the relative abundance of ParA and ParP. In contrast to figure 3, for the majority of the parameter space considered, there is only one biologically feasible steady state (which is mainly unstable). This region of the parameter space is highlighted in figure 9. The region of the parameter space where there exists three steady-state solutions has a low total amount of ParP compared with the total amount of ParA. The stability of the third steady state is always stable, whereas the second steady state is mainly unstable and the first steady state shows a stripe pattern which alters between stable and unstable. We highlight the fact that for the case in which both *A _{T}* and

*P*are very small, there exists only one steady state suggesting that in the protein conservation model, higher protein abundance creates additional steady states.

_{T}### 3.4. Numerical simulations of deterministic and stochastic models with cytoplasmic dimerization

We use the same methods of numerical simulation and initial conditions as we did for the protein conservation model without cytoplasmic dimerization. The results of the numerical simulations are presented in figure 10. The first row of figure 10 shows how the numerical solution of the protein conservation model with cytoplasmic dimerization evolves for both the deterministic and the stochastic cases in the low-protein abundance parameter regime (the total number of proteins is 13). Unlike the constant cytoplasmic pool model, we find the deterministic and stochastic trajectories converge to a single steady state (consistent with the findings of figure 9). We note that the stochastic model appears to converge to a steady state that has approximately equal abundances of ParP and total ParA—as was true for the constant cytoplasmic pool model and protein conservation model without cytoplasmic dimerization. The second row of figure 10 shows the numerical results for the medium protein abundance parameter regime (total protein number is 195). We find that in this case, both the constant cytoplasmic pool (QSS) model and the protein conservation model are in excellent agreement. Switching behaviour is found in the stochastic model and the deterministic model shows clear bistable behaviour. This model behaviour is very similar to the constant cytoplasmic pool model but was not found at all in the protein conservation model without cytoplasmic dimerization. Finally, the third row reveals the numerical solutions for the high-protein abundance regime (total number of proteins is 1950). For the most part, there is excellent agreement between the stochastic and the deterministic models. However, for one trajectory (produced using initial condition given by Goehring *et al*. [28]), the stochastic model converges to the high total ParA, low ParP steady state, whereas the deterministic model converges to the low total ParA, high ParP steady state. It appears that the noise levels were sufficiently high to kick the solution into another steady state—a phenomenon we did not observe in figure 5, which showed the results for the constant cytoplasmic pool model. By performing repeated simulations, we were able to find other examples of this behaviour in both the constant cytoplasmic pool model and the protein conservation model. These findings serve to show that the stochastic model solutions can differ from the corresponding deterministic solutions even in the high-protein abundance parameter regime. We discuss the implications of this in Discussion.

## 4. Discussion and future directions

Polarization is a fundamental biological process, whereby a cell specifies a spatial axis by localizing determinants such as proteins and lipids to opposite ends of the cell. Cell polarization is required for a variety of important cell functions including asymmetric division, directed migration and tissue integrity. Depending on the biological context, it may (or may not) be necessary for polarized cells to remain sensitive to additional signals. For instance, embryonic cells that rely on polarization for asymmetric division and cell fate specification will be relatively insensitive to further signals, which we refer to as hard polarization, while neutrophils that chase down invading bacteria need to quickly change their axis of polarity in response to chemical cues, which we term soft polarization. A potential regulatory mechanism used by cells to elicit different behaviours while using the same biochemical signalling pathways relies on tight control of active protein amount, either through transcriptional or biochemical regulation [29–31]. Indeed, different levels of protein expression can be used to differentiate between closely related cell types [32] and can be used as a diagnostic tool for several disease processes such as cancer and lupus [33–35].

In this paper, we demonstrate the influence of protein copy numbers on deterministic and stochastic models of Par protein asymmetry. To begin, we constructed a stochastic analogue of the Par protein model presented in R1 and compared it with its deterministic counterpart. Throughout the paper, we have referred to these models as ‘QSS’ or ‘constant cytoplasmic pool’ models, as the key assumption underlying these models is that the abundance of protein (be it anterior or posterior) in the cytoplasm is constant. To study the influence of protein abundance on these models, we varied the amount of protein in the cytoplasmic pools.

We explored the deterministic implementation of the constant cytoplasmic pool in detail and found that for the majority of the (*A _{y}*,

*P*) parameter space, three steady states existed, excluding regions where or . The dynamics of both the deterministic model and stochastic model are shown in figure 11. In other words, if the size of either cytoplasmic pool (anterior or posterior Par proteins) greatly exceeded the other, then the system exhibits only one steady state. Through a linear stability analysis of the steady-state equations and visualization in the phase plane portraits, we found that for regions of the parameter space producing three steady states, one steady state was always unstable while the other two were stable. These are characteristic properties of a bistable system. In the corresponding stochastic model, we explored three different ‘snapshots’ of the parameter space pertaining to ‘low’, ‘medium’ and ‘high’ protein abundances (table 1). Interestingly, we found that the cumulative mean of the stochastic model solution did not converge to the steady-state solution of the corresponding deterministic model for 'low’ and 'medium’ protein abundance regimes. For the ‘low’ protein abundance regime, we found that intrinsic noise dominated the stochastic solution and that regardless of initial conditions, the solution converged to approximately the same value (which had a non-zero amount of Total ParA and ParP), whereas the deterministic model exhibited bistable dynamics for the same parameters. In the ‘medium’ copy number parameter regime, the stochastic model exhibited switching behaviour in which the model trajectories jumped between the corresponding steady states of the deterministic model. In other words, the stochastic model solution was very sensitive to change where the solution appeared to converge to the corresponding deterministic steady state, the solution could then spontaneously jump to the other steady state. We hypothesize that in this protein abundance regime, a cell would be sensitive to external cues, with the spatial arrangement of cortex-bound Par proteins remaining very susceptible to change (what we refer to as ‘soft polarization’). This variation of behaviour between stochastic and deterministic models has been observed and studied in other biological systems such as forest insect populations [36] and neural firing rates [37]. Finally, we found that the stochastic and deterministic models were in close agreement for large protein abundances. We describe cells with high abundances of Par proteins as those which would undergo ‘hard polarization’. That is, the spatial arrangement of proteins would not be responsive to external cues.

_{y}To further study the role of protein abundance, we proposed a second model (considering both deterministic and stochastic versions) which constrained the total number of Par proteins. The constant cytoplasmic pools were replaced with dynamical pools which changed in size depending on whether proteins attached or detached from the cortex. By specifying initial conditions, we knew precisely how many anterior and posterior Par proteins were in our model ‘cell’ at any point in time. Interestingly, we found that without the inclusion of an additional dimerization reaction occurring in the cytoplasmic pool, bistable behaviour could not be observed for our protein conservation model (only a single steady state exists for both stochastic and deterministic models). However, we found that including such a dimerization reaction resulted in the appearance of bistable behaviour in both the deterministic and stochastic cases. The dynamics of the deterministic and stochastic versions of the model with cytoplasmic dimerization are shown in figure 12. We explored the deterministic model in detail, numerically determining the number of steady states for a large range of total ParA and ParP abundances and verifying the steady-state behaviours using linear stability analysis. The corresponding stochastic model was explored for snapshots of the protein abundance spectrum. With respect to the deterministic protein conservation model (including a cytoplasmic dimerization reaction), we found that bistable behaviour occurred when the total amount of anterior Par proteins far exceeded the total amount of posterior Par proteins. Unlike the constant cytoplasmic pool model, for low-protein abundances (13 Par proteins in total), the stochastic and deterministic models produce results in qualitative agreement (both converge to a single steady state) but not quantitative agreement (the stochastic model converges to a steady state with equal amounts of total ParA and ParP, while the deterministic model converges to a steady state with high ParP levels and zero total ParA levels). For medium protein abundances (195 Par proteins in total), stochastic switching behaviour was observed similar to that found in the constant cytoplasmic pool model while the deterministic implementation displayed bistability. Finally, in the case where there were 1950 proteins in total, we found reasonably good agreement between the stochastic and deterministic solutions. However, for one initial condition (close to the separatrix), the noise of the stochastic solution was enough to ‘kick’ the model to a steady state that was different from that of the corresponding deterministic model. This behaviour was exhibited in both the constant cytoplasmic pool and protein conservation models. Consistent with our notion of ‘hard polarization’, once a steady state was reached in the high-protein abundance regime, we found no evidence of any switching behaviour.

In the ‘medium’ copy number regime for both the QSS and conservation models, we found that trajectories can either behave like their mean-field counterparts and converge towards the corresponding steady state, or conversely, jump across the separatrix and converge towards the other mean-field basin of attraction. In other words, if noise is taken into account, then the separatrix becomes a barrier that the system may cross with finite probability. A well-studied theoretical bistable system which possesses this property is the Schloögl model and we refer the reader to [38] for further information. We note that, in general, switching behaviour in bistable systems can also be modelled as an activated process and, consequently, the corresponding transition rates can be modelled as functions of the Arrhenius type. In particular, there exists extensive literature [26,39] showing that, in the limit of low noise intensity, *σ*, the process of switching between two stable states is an activated process, i.e. the transition rate from one steady state to the other, *W*_{T}, is that such that *W*_{T} ∼ e^{−H} for , where *H* is a function of the parameters of the system and noise intensity. The consequence is that for low noise intensity, the probability of switching is exponentially suppressed.

Our model results can be tested via a series of experiments. By titrating the level of protein activity or expression, for instance using optogenetic controls [40], the response of a cell to noisy external signals can be quantified. Our model predicts that cells, such as motile cells, that generally tend to be sensitive to external signals would be less responsive to the same strength of external signal when polarity determinants are expressed at a higher amounts. Conversely, the model predicts that cells that are relatively insensitive to additional signals, such as embryonic cells, will more likely switch their polarity orientation if the amount of polarity determinants is reduced.

In this paper, we used a modelling approach to elucidate potential mechanisms regulating polarization behaviour in a cell, and in particular how a cell can adjust its sensitivity to a noisy environment by regulating the amount of relevant proteins. We used two distinct but complementary techniques: a stochastic modelling approach is suitable for situations where protein abundance is low, while a deterministic approach is suitable for high protein abundance. The two modelling approaches agree for high-copy numbers as would be expected according to the Law of Large Numbers. At low-copy number, our stochastic model is able to produce novel behaviours not possible with the deterministic model. Future work will extend this model comparison framework to yield a better understanding of the impact of noise, both internally and externally generated, on the ability of the model cell to assume a polarized state. To this extent, we will investigate probability density functions under a range of conditions. We present samples of such probability density functions in figures 13 and 14 for both the constant cytoplasmic pool model and conserved protein model with cytoplasmic dimerization in the medium copy number regime, due to the switching behaviour observed in the stochastic models. It is clear from both these figures that the models produce bimodal probability density functions. In the case of the constant cytoplasmic pool model, the non-zero states occur with far less probability than the zero ones—indicating that the cortex-bound species only occupy the cortex for brief periods. It is also worth noting that the ParP protein occupies the cortex with higher probability than the anterior Par molecules. The conserved protein model shows similar behaviour with the cytoplasmic species occupying non-zero states with higher probability. The sharp peaks in the cytoplasmic species reflect the high probability of the cortical-bound species being found in the zero state. These observations are consistent with the simulation results, indicating that the probability density functions are able to accurately capture aggregate behaviours of the stochastic models. An interesting multi-modal distribution appears for the cytoplasmic anterior Par protein species (*A _{y}*) which is due to the dimerization reaction. We reserve further discussion and analysis for a subsequent study.

We are also extending this model to a spatial setting, incorporating spatial movement of Par proteins on the cortex through diffusion to see if we can find the establishment of complementary domains with an analogous stochastic model. This model will account for the non-uniformly crowded cellular interior which is known to have an impact on both the diffusion coefficient and reaction rates [41–44]. We particularly wish to use this spatial model to further investigate our preliminary observations that cells that are likely to take on a soft polarization state have two-dimensional-type geometries where exchange between cytoplasmic and cortical compartments is negligible but where protein conservation plays a critical role [28,45], while hard polarization tends to occur in cells with distinct three-dimensional geometries where cytoplasmic diffusion, away from the cortex, may play a significant role such as in columnar epithelial cells or embryonic cells [6,46]. This leads our research into intriguing areas, asking whether cellular geometry may be playing a critical, but poorly understood, role in determining a polarized cell's response to additional signals.

## Funding statement

M.S. would like to thank the support from the Mathematical Biosciences Institute at The Ohio State University and NSF grant no. DMS0931642. A.T.D. is supported by NSF/NIH grant no. DMS/NIGMS-1361251.

- Received February 18, 2015.
- Accepted March 23, 2015.

- © 2015 The Author(s) Published by the Royal Society. All rights reserved.