**CDC-42 Interactions with Par Proteins Are Critical for Proper Patterning in Polarization**

**Sungrim Seirin-Lee 1,\*, Eamonn A. Gaffney <sup>2</sup> and Adriana T. Dawes 3,4,\***


Received: 10 August 2020; Accepted: 2 September 2020; Published: 5 September 2020

**Abstract:** Many cells rearrange proteins and other components into spatially distinct domains in a process called polarization. This asymmetric patterning is required for a number of biological processes including asymmetric division, cell migration, and embryonic development. Proteins involved in polarization are highly conserved and include members of the Par and Rho protein families. Despite the importance of these proteins in polarization, it is not yet known how they interact and regulate each other to produce the protein localization patterns associated with polarization. In this study, we develop and analyse a biologically based mathematical model of polarization that incorporates interactions between Par and Rho proteins that are consistent with experimental observations of CDC-42. Using minimal network and eFAST sensitivity analyses, we demonstrate that CDC-42 is predicted to reinforce maintenance of anterior PAR protein polarity which in turn feedbacks to maintain CDC-42 polarization, as well as supporting posterior PAR protein polarization maintenance. The mechanisms for polarity maintenance identified by these methods are not sufficient for the generation of polarization in the absence of cortical flow. Additional inhibitory interactions mediated by the posterior Par proteins are predicted to play a role in the generation of Par protein polarity. More generally, these results provide new insights into the role of CDC-42 in polarization and the mutual regulation of key polarity determinants, in addition to providing a foundation for further investigations.

**Keywords:** intracellular polarization; partial differential equations; sensitivity analysis

#### **1. Introduction**

Intracellular polarization, whereby a cell establishes a pattern and specifies a spatial axis by segregating proteins and other factors to distinct domains, is a fundamental and ubiquitous process. Polarization is implicated in a wide variety of biological phenomena including asymmetric cell division, cell migration, wound healing, and embryonic development [1–3]. Aberrant polarization is also thought to play a role in disease progression: a hallmark of the epithelial to mesenchymal transition (EMT) in malignant cells is the acquisition of a polarized migratory phenotype [4,5]. The same key polarity determinants, including the Par and Rho protein families, are required for forming the pattern associated with polarization in virtually all cell types and organisms [3]. Despite the importance of polarization and the highly conserved nature of the proteins involved, the mechanisms and signalling networks regulating this patterning process are not completely understood.

Two main protein families have been intensively studied for their role in polarization: Rho proteins and Par proteins (Figure 1). Rho proteins, also referred to as Rho GTPases, are monomeric G proteins that interconvert between an active GTP-bound state, and an inactive GDP-bound state. The active form of the protein can associate with the membrane at the periphery of the cell, while the inactive form is found diffusing in the cytoplasm [6]. Three members of this family, Rac, Rho and Cdc-42, have been extensively studied for their role in cell migration [6–9]. Polarized migrating cells establish front and back domains, with Cdc-42 and Rac segregating to the front, and Rho associating with the rear of the cell. Par proteins were first identified in the early embryos of the nematode worm *C. elegans* for their role in properly patterning the embryo prior to first division [3,10,11]. The asymmetric localization of the Par proteins specifies the anterior/posterior axis of the developing embryo, and is required for the asymmetric first division. PAR-3, PAR-6 and the atypical protein kinase aPKC bind to the membrane on the periphery of the cell and specify the anterior half of the cell, while PAR-1 and PAR-2 bind to the membrane on the posterior half. Loss of polarity in the first cell cycle is lethal to the embryo [12]. Recent experimental observations suggest that Par and Rho proteins rely on mutual interaction and feedback to produce the pattern of proteins associated with polarization: T cells require both Cdc-42 and Par proteins for polarization [13], and *C. elegans* embryos require CDC-42 for proper patterning of PAR-2 and PAR-6 [14,15]. Despite their involvement in polarization, the dynamics of Par and Rho proteins have been studied largely independently from each other.

**Figure 1.** Interaction network of PAR proteins and the Rho GTPase CDC-42. (**A**) Top: Interaction network between PAR proteins and CDC-42 consistent with experimental results. References for interactions (i)–(viii) are given in Section 2.1 of the main text. Bottom: Interaction network used in this study, with cytosol to membrane translocation represented by activation, and membrane to cytosol translocation represented by inhibition. (**B**) Reduction of interaction network in (**A**) into subnetworks. The interaction network is separated into three modules: (b1) shows the interactions between the anterior Par proteins PAR-6 and PKC-3 and the posterior Par proteins PAR-1 and PAR-2, (b2) shows the interactions between the anterior Par protein PAR-3 and the posterior Par proteins PAR-1 and PAR-2, and (b4) shows the interactions between the anterior and posterior Par proteins and the Rho GTPase CDC-42. Note that the upregulation of PAR-6 by PAR-3 is implicitly captured by the inhibitory action of PAR-3 on PAR-2, since the latter in turn downregulates PAR-6. As described in the text, (b1) and (b2) can be reduced to (b3). The green dotted arrows indicates two sequential inhibition interactions (⊥) that are equivalent to an activating interaction (↑, solid green line). (**C**) Network (b3) merged with (b4) is the fundamental network investigated here. Rate constants *q<sup>a</sup>* <sup>3</sup>, *<sup>q</sup><sup>c</sup>* <sup>3</sup>, *k p* <sup>3</sup> , *<sup>k</sup><sup>a</sup>* <sup>3</sup>, and *<sup>k</sup><sup>c</sup>* <sup>3</sup> appear in the corresponding model Equations (4).

While polarization is observed in many different systems, we focus on protein patterning in the early *C. elegans* embryo [16] (Figure 2). The establishment phase, which is initiated by fertilization, moves symmetrically distributed Par proteins, including PAR-3, PAR-6 and aPKC, to the anterior half of the cell. The cleared area in the posterior half is then available for binding by the posterior Par proteins, including PAR-1 and PAR-2. Once the polarized domains are established, the cell transitions into the maintenance phase, where it maintains the asymmetric protein pattern as the cell prepares for first division. While the protein transport associated with polarization is important for establishment of polarization [17], we focus on the role of biochemical interactions in the generation and maintenance of polarized domains, independent of advective flow. Members of the Rho protein family, notably CDC-42, are thought to be important for polarization in the early embryo, but how they interact with and regulate Par proteins is not clear.

**Figure 2.** Schematic of model components and geometry. (**A**) Model reduction from three dimensions to one dimension in space. A cross section of the bulk cytosol region is coloured by red, the surface is coloured by gray, with *uc*, *um* denoting generic cytosolic and membrane densities, respectively. Fast cytosolic diffusion timescales compared to the timescales of interest lead to effectively homogeneous cytosolic concentrations, that can be captured by membrane concentrations via conservation principles, as discussed in Section 2.2. This allows the model to be reduced to membrane equations, with axisymmetry about the major axis of the cell for a simple geometry yielding the final simplification, as detailed in Appendix A.1. (**B**) Representative simulations for wild type case and *cdc-42(RNAi)* case. Simulations start with initial conditions given by initial conditions, Equation (5). After 30 min, the solutions shown have reached steady state. Here, and throughout the paper, concentrations are in units of M/*L*, where M is the characteristic protein number of Table A1, and *L* is the perimeter of the projected membrane in (**A**). (**C**) The range of values for the parameters *k p* <sup>3</sup> , and *<sup>k</sup><sup>a</sup>* <sup>3</sup> that satisfy both wild type and *cdc-42(RNAi)* behaviours which serves as a validation of the model against observations. The red circles indicate the parameter region matching both wild type and *cdc-42(RNAi)* behaviours. The black square indicates the range of (*k p* <sup>3</sup> , *<sup>k</sup><sup>a</sup>* <sup>3</sup>) used for sensitivity analysis.

In this paper, we develop a continuum model of Par and Rho protein dynamics in the generation and maintenance of polarization in the early *C. elegans* embryo from the rational simplification for the underlying interactions among CDC-42 and the PAR proteins that have been uncovered by numerous experimental studies. This model also utilizes the large ratio of diffusive transport scales between the cell cytosol and membrane and a simple representation of the cell geometry. We demonstrate that the resulting partial differential equation model is consistent with observations of CDC-42, highlighting the requirement for CDC-42 during the maintenance of polarization in the early *C. elegans* [18,19]. To elucidate the detailed mechanism of interactions capable of reproducing experimental observations, we perform sensitivity analysis and minimal network analysis to identify and characterize predictions for key cross talk and mutual regulation interactions among CDC-42 and the PAR proteins in controlling the generation and maintenance of cellular polarization. Due to the

conserved nature of both the process and proteins involved in polarization, the insights gained in this study apply broadly to other biological systems beyond *C. elegans*.

#### **2. The Mathematical Model**

#### *2.1. Network of Par Protein and Cdc-42 Interactions*

To investigate the role potential role of interactions between Par and Rho protein family members in the generation and maintenance of distinct spatial domains, we construct a network of interactions that is consistent with experimental data. Justification for each arrow in Figure 1A is as follows:


Figure 1A shows a full schematic diagram and the network of interactions corresponding to the descriptions (i)–(viii). Some of these interactions may not be direct, but are mediated by other proteins. For instance, enhanced cortical dissociation of the anterior Par proteins by PAR-2 is likely mediated by PAR-1 [27]. In addition, the role of CDC-42 and its interactions with the Par proteins are not entirely clear, and some observations are not accounted for in our model, such as the appearance of CDC-42 on the cortex even in the absence of PAR-6 [14]. Nevertheless, experimental observations of all interactions listed above have been reported in the literature.

To determine the fundamental network of interactions based on the above experimental evidence, we separated the interaction network in Figure 1A into three subnetworks, Figure 1B(b1,b2,b4). In Figure 1B(b1), we consider the interactions between [PAR-6, PAR-2 and PAR-1] and [PKC-3, PAR-2 and PAR-1]. In this subnetwork, the inhibition of PAR-1 by PAR-6 and inhibition of PAR-6 by PAR-2 (dotted green line) is equivalent to PAR-2 activation of PAR-1 (solid green line). Since PAR-6 and PKC-3 only dissociate PAR-1 and PAR-2 when PAR-6 and PKC-3 are bound together in a complex, we can reduce this subnetwork to a mutual inhibition between [PAR-6 and PKC-3] and [PAR-1 and PAR-2] (Figure 1B(b1)). In a similar manner, the subnetwork containing [PAR-3, PAR-2 and PAR-1] can be reduced to a network of mutual inhibition between [PAR-3] and [PAR-1 and PAR-2] (Figure 1B(b2)). We further simplify these subnetworks by grouping the Par proteins according to their localization: anterior Par proteins (aPAR) and posterior Par proteins (pPAR) (Figure 1B(b3)). Finally, we combine the mutual Par protein inhibition subnetwork with the [CDC-42, PAR-6/PKC-3 and PAR-2] subnetwork (Figure 1B(b4)) to arrive at the fundamental interaction network of [aPAR, pPAR, and CDC-42] (Figure 1C).

#### *2.2. Model Equations for Apar, Ppar and Cdc-42 Network*

We define the cytosol by <sup>Ω</sup> <sup>⊂</sup> <sup>R</sup><sup>3</sup> and the membrane by *<sup>∂</sup>*Ω(<sup>≡</sup> <sup>Γ</sup>) so that <sup>Ω</sup> <sup>∪</sup> *<sup>∂</sup>*<sup>Ω</sup> represents *C. elegans* single cell embryo (Figure 2A). We represent the concentrations of anterior proteins (aPARs, i.e., PAR-6, PAR-3, and PKC-3) in the membrane and the cytosol by [*Am*(**x**, *t*)] and [*Ac*(**x**, *t*)], respectively, and the concentrations of posterior proteins (pPARs, i.e., PAR-1 and PAR-2) in the membrane and the cytosol by [*Pm*(**x**, *t*)] and [*Pc*(**x**, *t*)], respectively, and the concentrations of CDC-42

in the membrane and the cytosol by [*Cm*(**x**, *<sup>t</sup>*)] and [*Cc*(**x**, *<sup>t</sup>*)], respectively, where **<sup>x</sup>** <sup>∈</sup> <sup>R</sup><sup>3</sup> and *<sup>t</sup>* <sup>∈</sup> [0, <sup>∞</sup>). The dynamics of the network shown in Figure 1C can then be described as follows:

$$\begin{aligned} \frac{\partial [A\_{\mathfrak{m}}]}{\partial t} &= D\_{\mathfrak{m}}^{A} \nabla\_{\Gamma}^{2} [A\_{\mathfrak{m}}] + \left\{ \gamma\_{\mathfrak{a}} + F\_{\mathrm{on}}^{A} ([\mathbf{C}\_{\mathfrak{m}}]) \right\} [A\_{\mathfrak{c}}] - \left\{ \mathbf{a}\_{\mathfrak{d}} + F\_{\mathrm{off}}^{A} ([P\_{\mathfrak{m}}]) \right\} [A\_{\mathfrak{m}}] & \quad \text{on } \mathbf{x} \in \partial \Omega, \\\ \partial [A\_{\mathfrak{c}}] &= D^{A} \nabla^{2} [A\_{\mathfrak{m}}] \end{aligned}$$

$$\begin{aligned} \frac{\partial \{\mathcal{A}\}}{\partial t} &= D\_{\varepsilon}^{A} \nabla^{2} [A\_{\varepsilon}] & \text{on } \mathbf{x} \in \Omega, \\ D\_{\varepsilon}^{A} \frac{\partial [A\_{\varepsilon}]}{\partial \mathbf{n}} &= -\{\gamma\_{a} + F\_{\text{on}}^{A}([\mathbb{C}\_{\mathcal{m}}])\}[A\_{\varepsilon}] + \{a\_{\mathfrak{q}} + F\_{\text{off}}^{A}([P\_{\mathcal{W}}])\}[A\_{\mathcal{W}}] & \text{on } \mathbf{x} \in \partial \Omega, \end{aligned}$$

$$\begin{aligned} \frac{\partial[P\_{\rm m}]}{\partial t} &= D\_{\rm m}^{P} \nabla\_{\Gamma}^{2}[P\_{\rm m}] + \gamma\_{p}[P\_{\rm c}] - \{a\_{\mathcal{P}} + F\_{\rm off}^{P}([A\_{\rm m}])\}[P\_{\rm m}] & & \text{on } \mathbf{x} \in \partial \Omega, \\ \frac{\partial[P\_{\rm c}]}{\partial \mathbf{x}} &= D\_{\rm c}^{P} \nabla^{2}[P\_{\rm m}] & & \text{on } \mathbf{x} \in \Omega, \end{aligned} \tag{1}$$

$$\begin{aligned} \frac{\partial^{\left[\mathbf{I}\cdot\mathbf{I}\right]}}{\partial t} &= D\_{\mathbf{c}}^{P} \nabla^{2} [P\_{\mathbf{m}}] & \text{on } \mathbf{x} \in \Omega, \\ D\_{\mathbf{c}}^{P} \frac{\partial [P\_{\mathbf{c}}]}{\partial \mathbf{n}} &= -\gamma\_{P} [P\_{\mathbf{c}}] + \left\{ a\_{\mathbf{d}} + F\_{\mathrm{off}}^{P} ([A\_{\mathbf{m}}]) \right\} [P\_{\mathbf{m}}] & \text{on } \mathbf{x} \in \partial \Omega, \end{aligned}$$

$$\begin{cases} \frac{\partial [\mathbb{C}\_{m}]}{\partial t} = D\_{m}^{\mathbb{C}} \nabla\_{\mathbb{T}}^{2} [\mathbb{C}\_{m}] + \gamma\_{\mathbb{c}} [\mathbb{C}\_{\mathbb{c}}] - \{a\_{\mathbb{c}} + F\_{\text{off}}^{\mathbb{C}}([P\_{m}], [A\_{m}])\} [\mathbb{C}\_{m}] & \text{on } \mathbf{x} \in \partial \Omega, \\\frac{\partial [\mathbb{C}\_{\mathbb{c}}]}{\partial t} = D\_{\mathbb{c}}^{\mathbb{C}} \nabla^{2} [\mathbb{C}\_{\mathbb{c}}] & \text{on } \mathbf{x} \in \Omega, \\\mathbb{C}^{\mathbb{N}} \mathbb{C}^{\mathbb{L}} & \end{cases}$$

$$D\_{\varepsilon}^{\mathbb{C}} \frac{\partial [\mathbb{C}\_{\varepsilon}]}{\partial \mathbf{n}} = -\gamma\_{\varepsilon} [\mathbb{C}\_{\varepsilon}] + \{a\_{\varepsilon} + F\_{\mathrm{off}}^{\mathbb{C}}([P\_{m}], [A\_{m}])\} [\mathbb{C}\_{m}] \tag{10.4} \qquad \text{on } \mathbf{x} \in \partial \Omega,$$

where *D<sup>A</sup> <sup>m</sup>*, *D<sup>P</sup> <sup>m</sup>* and *D<sup>C</sup> <sup>m</sup>* are the diffusion coefficients of aPARs, pPARs and CDC-42 in the membrane, respectively, and *D<sup>A</sup> <sup>c</sup>* , *D<sup>P</sup> <sup>c</sup>* and *D<sup>C</sup> <sup>c</sup>* are the diffusion coefficients of aPARs, pPAR and CDC-42 in the cytosol, respectively. **n** is the inner normal vector on *∂*Ω. We assume that *F<sup>A</sup>* on, *F<sup>A</sup>* off and *<sup>F</sup><sup>P</sup>* off are increasing functions of CDC-42 concentration, pPAR concentration, and aPAR concentration, respectively. We also assume that *F<sup>C</sup>* off is an increasing function of pPAR concentration but a decreasing function of aPAR concentration. For the purposes of this study, we choose the following functional forms [22,28];

$$\begin{split} F\_{\text{on}}^{A}([\mathbb{C}\_{\text{m}}]) &= \frac{q\_{3}^{a}[\mathbb{C}\_{\text{m}}]^{2}}{q\_{1}^{a} + q\_{2}^{a}[\mathbb{C}\_{\text{m}}]^{2}}, \quad F\_{\text{off}}^{A}([P\_{\text{m}}]) = \frac{k\_{3}^{a}[P\_{\text{m}}]^{2}}{k\_{1}^{a} + k\_{2}^{a}[P\_{\text{m}}]^{2}},\\ F\_{\text{off}}^{P}([A\_{\text{m}}]) &= \frac{k\_{3}^{p}[A\_{\text{m}}]^{2}}{k\_{1}^{p} + k\_{2}^{p}[A\_{\text{m}}]^{2}}, \quad F\_{\text{off}}^{C}([P\_{\text{m}}]\_{\text{\textprime}}[A\_{\text{m}}]) = \frac{k\_{3}^{c}[P\_{\text{m}}]^{2}}{k\_{1}^{c} + k\_{2}^{c}[P\_{\text{m}}]^{2}} + \frac{q\_{3}^{c}}{q\_{1}^{c} + q\_{2}^{c}[A\_{\text{m}}]^{2}}. \end{split} \tag{2}$$

Rate constants of the form *k j <sup>i</sup>* denote interactions that include pPAR while rate constants of the form *q j i* denote interactions strictly between aPAR and CDC-42.

Parameter values corresponding to the final simplified model are summarized in Table A1. Note that some parameters and variables are redefined during this simplification process. The first step in this simplification exploits differences in diffusion dynamics. We require *D<sup>A</sup> <sup>c</sup> <sup>D</sup><sup>A</sup> <sup>m</sup>*, *D<sup>P</sup> <sup>c</sup> <sup>D</sup><sup>P</sup> <sup>m</sup>* and *D<sup>C</sup> <sup>c</sup> <sup>D</sup><sup>C</sup> <sup>m</sup>* since diffusion in the cytosol is much faster than diffusion in the membrane [29,30]. Furthermore, this fast cytosolic diffusion, if sufficiently large, results in homogeneous spatial distribution of cytoplasmic proteins over the time scale of interest (Figure A1) and allows us to reduce the model (1) to a shadow system. By taking *D<sup>A</sup> <sup>c</sup>* , *D<sup>P</sup> <sup>c</sup>* , *D<sup>C</sup> <sup>c</sup>* → ∞, the leading order approximations for cytosolic protein concentrations quickly approach a homogeneous steady states such that

$$A\_{\varepsilon}(\mathbf{x},t) \to \frac{1}{|\Omega|} \int\_{\Omega} [A\_{\varepsilon}](\mathbf{x},t) d\mathbf{x}, \quad P\_{\varepsilon}(\mathbf{x},t) \to \frac{1}{|\Omega|} \int\_{\Omega} [P\_{\varepsilon}](\mathbf{x},t) d\mathbf{x}, \quad \mathbb{C}\_{\varepsilon}(\mathbf{x},t) \to \frac{1}{|\Omega|} \int\_{\Omega} [\mathbb{C}\_{\varepsilon}](\mathbf{x},t) d\mathbf{x}.$$

The following conservation relations hold:

$$\begin{aligned} \frac{d}{dt} \left( \int\_{\Omega} [A\_c](\mathbf{x}, t) d\mathbf{x} + \int\_{\Gamma} [A\_m](\mathbf{x}, t) d\mathbf{x} \right) &= 0, \\ \frac{d}{dt} \left( \int\_{\Omega} [P\_c](\mathbf{x}, t) d\mathbf{x} + \int\_{\Gamma} [P\_m](\mathbf{x}, t) d\mathbf{x} \right) &= 0, \\ \frac{d}{dt} \left( \int\_{\Omega} [\mathbf{C}\_c](\mathbf{x}, t) d\mathbf{x} + \int\_{\Gamma} [\mathbf{C}\_m](\mathbf{x}, t) d\mathbf{x} \right) &= 0, \end{aligned}$$

and the total mass of the model (1) is conserved. We define the conserved total concentrations as follows:

$$\begin{aligned} A\_{\text{tot}} & \equiv \int\_{\Omega} [A\_{\text{c}}](\mathbf{x}, t) d\mathbf{x} + \int\_{\Gamma} [A\_{\text{m}}](\mathbf{x}, t) d\mathbf{x}, \qquad P\_{\text{tot}} \equiv \int\_{\Omega} [P\_{\text{c}}](\mathbf{x}, t) d\mathbf{x} + \int\_{\Gamma} [P\_{\text{m}}](\mathbf{x}, t) d\mathbf{x}, \\ \mathsf{C}\_{\text{tot}} & \equiv \int\_{\Omega} [\mathsf{C}\_{\text{c}}](\mathbf{x}, t) d\mathbf{x} + \int\_{\Gamma} [\mathsf{C}\_{\text{m}}](\mathbf{x}, t) d\mathbf{x}, \end{aligned}$$

giving us

$$\begin{split} A\_{\varepsilon}(\mathbf{x},t) &\approx \frac{1}{|\Omega|} \left( A\_{tot} - \int\_{\Gamma} [A\_{m}](\mathbf{x},t) d\mathbf{x} \right), \quad P\_{\varepsilon}(\mathbf{x},t) \approx \frac{1}{|\Omega|} \left( P\_{tot} - \int\_{\Gamma} [P\_{m}](\mathbf{x},t) d\mathbf{x} \right), \\ C\_{\varepsilon}(\mathbf{x},t) &\approx \frac{1}{|\Omega|} \left( \mathbb{C}\_{tot} - \int\_{\Gamma} [\mathbb{C}\_{m}](\mathbf{x},t) d\mathbf{x} \right), \end{split}$$

at leading order in an asymptotic approximation based on the cytosolic diffusion dominating all other possible diffusive scales in the model. Thus, working to this level of approximation, model (1) is reduced to the surface model:

$$\begin{cases} \frac{\partial [A\_m]}{\partial t} = D\_m^A \nabla\_\Gamma^2 [A\_m] + \{\gamma\_\ell + F\_{\text{off}}^A([\mathcal{C}\_m])\} \left( \frac{A\_{\ell\ell t}}{|\Omega|} - \frac{1}{|\Omega|} \int\_\Gamma [A\_m] d\mathbf{x} \right) - \{a\_d + F\_{\text{off}}^A([P\_m])\} [A\_m] & \text{on } \mathbf{x} \in \Gamma, \\\frac{\partial [P\_m]}{\partial t} = D\_m^P \nabla\_\Gamma^2 [P\_m] + \gamma\_\ell \left( \frac{P\_{\text{flat}}}{|\Omega|} - \frac{1}{|\Omega|} \int\_\Gamma [P\_m] d\mathbf{x} \right) - \{a\_p + F\_{\text{off}}^P([A\_m])\} [P\_m] & \text{on } \mathbf{x} \in \Gamma, \\\frac{\partial [\mathcal{C}\_m]}{\partial t} = D\_m^C \nabla\_\Gamma^2 [\mathcal{C}\_m] + \gamma\_\ell \left( \frac{\mathcal{C}\_{\ell\ell}}{|\Omega|} - \frac{1}{|\Omega|} \int\_\Gamma [\mathcal{C}\_m] d\mathbf{x} \right) - \{a\_\ell + F\_{\text{off}}^C([P\_m], [A\_m])\} [\mathcal{C}\_m] & \text{on } \mathbf{x} \in \Gamma. \end{cases} \tag{3}$$

We define the model geometry as the surface of a solid of revolution found by rotating the arc of the membrane about the *AP*-axis (Figure 2A). Since we are only interested in the dynamics on the membrane, we seek to reduce model (3) to a one dimensional model on *x* ∈ [0, *L*] where *x* is the arclength along the perimeter of the projected membrane, depicted in Figure 2A, *x* = 0 corresponds to the leftmost point on the membrane and *L* is the cell perimeter. This reduction to a 1D domain also ensures the minimal network and eFAST analysis can be accomplished in a reasonable amount of time. Periodic boundary conditions are imposed at *x* = 0, *L* (Figure 2A). We further approximate the geometry effectively by a cylinder of radius *H L* with asymptotically short caps and we define the tilde variables and parameters via

$$\bar{A}\_{\mathfrak{m}} = \pi H A\_{\mathfrak{m}}, \ \bar{P}\_{\mathfrak{m}} = \pi H P\_{\mathfrak{m}}, \ \bar{\mathsf{C}}\_{\mathfrak{m}} = \pi H \mathsf{C}\_{\mathfrak{m}}, \ \bar{\gamma}\_{\mathfrak{a}}/\gamma\_{\mathfrak{a}} = \bar{\gamma}\_{\mathfrak{p}}/\gamma\_{\mathfrak{p}} = \bar{\gamma}\_{\mathfrak{c}}/\gamma\_{\mathfrak{c}} = \bar{q}\_{3}^{\mathfrak{a}}/q\_{3}^{\mathfrak{a}} = \mathrm{L}\pi H / |\Omega|.$$

As detailed in Appendix A.1, we rewrite Equation (3) in terms of these new parameters, together with the assumption of axisymmetry and geometrical approximations that are stated and justified in Appendix A.1. On dropping tildes, we have

$$\begin{split} \frac{\partial[A\_{m}]}{\partial t} &= D\_{m}^{A} \frac{\partial^{2}}{\partial \mathbf{x}^{2}} [A\_{m}] + \left\{ \gamma\_{a} + F\_{\text{on}}^{A}([\mathbb{C}\_{m}]) \right\} \left( \frac{A\_{\text{tot}}}{L} - \frac{1}{L} \int\_{0}^{L} [A\_{m}] dx \right) - \left\{ a\_{\text{t}} + F\_{\text{off}}^{A}([P\_{\text{m}}]) \right\} [A\_{\text{m}}],\\ \frac{\partial[P\_{m}]}{\partial t} &= D\_{m}^{P} \frac{\partial^{2}}{\partial \mathbf{x}^{2}} [P\_{m}] + \gamma\_{p} \left( \frac{P\_{\text{tot}}}{L} - \frac{1}{L} \int\_{0}^{L} [P\_{m}] dx \right) - \left\{ a\_{p} + F\_{\text{off}}^{P}([A\_{m}]) \right\} [P\_{m}],\\ \frac{\partial[\mathbb{C}\_{m}]}{\partial t} &= D\_{m}^{\mathbb{C}} \frac{\partial^{2}}{\partial \mathbf{x}^{2}} [\mathbb{C}\_{m}] + \gamma\_{c} \left( \frac{\mathbb{C}\_{tot}}{L} - \frac{1}{L} \int\_{0}^{L} [\mathbb{C}\_{m}] dx \right) - \left\{ a\_{c} + F\_{\text{off}}^{C}([P\_{m}], [A\_{m}]) \right\} [\mathbb{C}\_{m}]. \end{split} \tag{4}$$

The parameter values in Table A1 are with respect to this final model (Equation (4)). The model was simulated by custom software written in C. The PDEs were solved via an implicit numerical scheme using standard finite-difference methods. Source code used to generate the results in this paper are available upon request.

#### *2.3. Parameter Values*

Before symmetry breaking in the *C. elegans* embryo, aPAR is spatially homogeneously distributed on the membrane and pPAR is spatially homogeneously distributed in the cytosol. Gotta et al [19] experimentally demonstrated that loss of CDC-42 by RNAi results in a loss of polarity, with low PAR-6 and high PAR-2 levels on the membrane. Thus, we choose representative kinetic parameters such that aPar, pPar and CDC-42 establish distinct spatial domains under wild type conditions, but fail to polarize when CDC-42 is absent (Figure 2B). For sensitivity analysis, we restrict our parameters as shown in Figure 2C, corresponding to parameter values that are consistent with both wild type and *cdc-42(RNAi)* experimental observations. The fact that such a large range of parameter values produces the appropriate model behaviours provides further validation of the model's predicted importance for CDC-42 in polarization. Parameter names and definitions for the final model, Equation (4), are summarized in Table A1, together with the non-dimensionalization used in the numerical investigation of the model.

#### *2.4. Initial Conditions*

We use two sets of initial conditions to evaluate either maintenance (Section 3.1.1) or generation (Section 3.2.1) of polarization. For simulating maintenance of polarization, we specify polarized initial conditions with aPAR and CDC-42 high in the anterior (left and right) region of the domain and low in the posterior (middle) region, and the reverse profile for pPAR (high in the posterior, low in the anterior). See Figure 2 for a schematic of the model geometry and the anterior/posterior regions. The high and low values are derived from the stationary long time asymptotic solution of Equation (4) (Figure 2B).

For simulating generation of polarization, we specify initial conditions with a small spatial perturbation:

$$\left[P\_{\mathfrak{m}}\right](\mathbf{x},0) = \epsilon\_p \delta(\mathbf{x} - \mathbf{L}/2), \left[A\_{\mathfrak{m}}\right](\mathbf{x},0) = \frac{A\_{\mathrm{tot}}}{L}(1 + \epsilon \Phi\_{\mathfrak{k}}(\mathbf{x})), \left[\mathbb{C}\_{\mathfrak{m}}\right](\mathbf{x},0) = \mathbb{C}\_{0}(1 + \epsilon \Phi\_{\mathfrak{c}}(\mathbf{x})), \tag{5}$$

where *x* ∈ [0, *L*], *δ* is the delta function, and *C*<sup>0</sup> is the equilibrium initial concentration of CDC-42. *φa*(*x*) and *φc*(*x*) are perturbation functions, *<sup>p</sup>* is the strength of the initial external perturbation signal and is the magnitude of the perturbation.

#### *2.5. Minimal Network Analysis*

To determine the minimal set of interactions required for maintenance (Section 3.1.1) or generation (Section 3.2.1) of a polarization pattern, we devised a method of minimal network analysis. Starting from the fundamental network containing all interactions (Figure 1C), we first remove individual interactions from the network. This reduced network is simulated and evaluated for the presence of a pattern using initial conditions described above. The total simulation time for pattern evaluation is 16 min for maintenance of polarization and 30 min for generation of polarization. While the numerical solution may not have converged to a stationary solution by the time of evaluation, it is clear whether a pattern has been initiated. We determine which interactions are common to all reduced networks that are capable of maintaining/generating polarization. A network consisting of these common interactions is tested for its capacity to maintain/generate a polarized pattern. If the common interaction network is not capable of maintaining/generating a pattern consistent with polarization, interactions are reintroduced into the model, first individually and then in pairs, and the model is simulated and evaluated for the presence of a polarized pattern. The reduced network(s) with the least number of interactions that is capable of maintaining/generating polarization is retained for further analysis.

#### *2.6. eFAST Sensitivity Analysis*

Sensitivity analysis is used to determine the influence of model parameters on the dynamics of the network via the use of Extended Fourier Amplitude Sensitivity Test (eFAST) [31,32]. In brief, we define sensitivity functions whose dynamics reveals self-organization and track their changes in response to changes in parameter values. In particular, the sensitivity functions for aPAR, pPAR, and CDC-42, denoted by *F<sup>A</sup> <sup>S</sup>* , *<sup>F</sup><sup>P</sup> <sup>S</sup>* and *<sup>F</sup><sup>C</sup> <sup>S</sup>* , respectively, are defined as

$$\begin{aligned} F\_S^A &= \sqrt{\int\_0^L \left| \frac{[A\_m](\mathbf{x}, t^\*) - [A\_m](\mathbf{x}, 0)}{t^\*} \right|^2}, \\ F\_S^P &= \sqrt{\int\_0^L \left| \frac{[P\_m](\mathbf{x}, t^\*) - [P\_m](\mathbf{x}, 0)}{t^\*} \right|^2}, \\ F\_S^C &= \sqrt{\int\_0^L \left| \frac{[\mathbf{C}\_m](\mathbf{x}, t^\*) - [\mathbf{C}\_m](\mathbf{x}, 0)}{t^\*} \right|^2} \end{aligned} \tag{6}$$

where *t* ∗ is the time scale of interest. This time scale is chosen as 16 min for maintenance of polarization, and 30 min for generation of polarization, consistent with experimental observations [33,34]. These sensitivity functions give a quantitative measure of how much the polarity pattern (the model output) has been altered in response to changes in parameter values (the model input). We calculate two sensitivity measures, the first order index (*Si*) and total-effect index (*STi* ). These measures are defined as follows:

$$S\_i \equiv \frac{\text{Variance of the expected model output } y \text{ with respect to parameter } p\_i}{\text{Total variance}} = \frac{V\_i}{V\_{total}},$$

$$S\_{T\_i} \equiv \text{Total effect (contribution) of parameter } p\_i \text{ to the output variance} = 1 - \frac{V\_{-i}}{V\_{total}}.$$

where *V*−*<sup>i</sup>* is the effect of any order that does not include the factor *i*. The first order index indicates the influence of parameter *pi* on the variance of the polarization measure (the model output), independent of interactions with the other parameters. The total-effect index indicates the effect of parameter *pi* when interactions with the other parameters are included. These two measures give a full quantification of the importance of parameter *pi* and whether the extent whether this is a direct influence or through interactions with other parameters or both. In this study, we focus on the parameters *q<sup>a</sup>* <sup>3</sup>, *<sup>q</sup><sup>c</sup>* <sup>3</sup>, *k p* <sup>3</sup> , *<sup>k</sup><sup>a</sup>* <sup>3</sup> and *<sup>k</sup><sup>c</sup>* <sup>3</sup> which determine the magnitude of the interaction functions, Equation (2).

A more in depth discussion of the eFAST sensitivity method and its detailed accompanying calculations can be found in Appendix A.3. The values used for sensitivity analysis are given in Table A2.

#### **3. Results**

#### *3.1. Critical Network Interactions and Parameters for Maintenance of Par Protein Polarization*

It has been observed that CDC-42 is required during maintenance but not establishment of polarization in the early *C. elegans* embryo [18,19]. However, it is not clear how CDC-42 is interacting with or regulating the anterior and posterior Par proteins to ensure maintenance of polarization in the early embryo.

#### 3.1.1. Minimal Network for Maintenance of Par Protein Polarization

Using our method of minimal network analysis, we aim to determine the minimal interaction network between anterior Par proteins (aPAR), CDC-42 and posterior Par proteins (pPAR) that may maintain spatial polarization. As shown in Figure 3A, the activation of aPAR by CDC-42 ((a1), *q<sup>a</sup>* <sup>3</sup> interaction) and inhibition of pPAR by aPAR ((a3), *k p* <sup>3</sup> interaction) are always present in reduced networks that are able to maintain polarization, suggesting these interactions are critical. The ability of the model to polarize in the absence of the other interactions ((a2), (a4), (a5)) show these interactions are less critical for polarity maintenance.

A minimal network consisting of only these two interactions (*q<sup>a</sup>* <sup>3</sup> and *k p* <sup>3</sup> , Figure 3C) is not sufficient to maintain polarization. As shown in Figure 3B, adding CDC-42 inhibition by pPAR ((b1), *k<sup>c</sup>* <sup>3</sup> interaction) or aPAR inhibition by pPAR ((b2), *<sup>k</sup><sup>a</sup>* <sup>3</sup> interaction) back to the network also cannot maintain polarization. Adding CDC-42 activation by aPAR ((b3), *q<sup>c</sup>* <sup>3</sup> interaction) allows polarity maintenance, although pPAR does not vary much over the domain. Further reinforcing the importance of CDC-42 for maintenance of polarity, we find that polarization of aPAR can be maintained with only two interactions (Figure 3D, *q<sup>a</sup>* <sup>3</sup> and *k p* <sup>3</sup> networks) and without mutual inhibition by pPAR. This is consistent with previous results, showing that aPARs intially polarize during the establishment phase, transiently establishing an anterior Par protein domain. However, this polarization can not be maintained, and the aPAR domain gradually creeps back toward the posterior pole, resulting in an eventual loss of polarity [17,18]. Our model reproduces the initial aPAR polarization in the absence of pPAR.

Taken together, the results suggest that the minimal network that is capable of maintaining polarization of aPAR and pPAR involves mutual activation of aPAR and CDC-42 and aPAR-mediated inhibition of pPAR (Figure 3B(b3)).

**Figure 3.** Minimal network needed for maintenance of Par protein polarization. (**A**–**D**) Representative simulations for the effects of network interactions. Dotted gray lines indicate interactions excluded from the network. Simulations are shown at *t* = 16 min, the approximate amount of time the early *C. elegans* embryo spends in the maintenance phase [33]. Initial conditions are given by the corresponding stationary solution (Figure 2B). (**A**) (a1)–(a5) Simulation results of networks omitting a single interaction. Minimal networks that produce the correct spatial pattern and are retained for further analysis are boxed in gray. (**B**) Red arrows indicate interactions identified in (**A**) that are retained in these additional networks. (b1)–(b3) Simulation results of networks omitting two interactions. (**E**,**F**) Variance-based sensitivity analysis results using eFAST. Parameters with higher values have a more significant effect on the network dynamics. (**E**) The first order index (the effect of a parameter on the model dynamics independent of the other parameters) and (**F**) The total effect index (the effect of a parameter including the interactions with other all model parameters) both indicate that CDC-42 interactions with aPAR are critical for producing the correct spatial pattern.

#### 3.1.2. Critical Parameters for Maintenance of Par Protein Polarization

To determine the effect of the network interactions on maintenance of Par protein polarity with the fundamental network, we performed a sensitivity analysis using the eFAST method [32,35]. See Appendix A.3 for more details about this approach. We wish to compare the critical interactions found by minimal network analysis with the important parameters indicated by the sensitivity analysis. We evaluated the first order index, *Si*, for the five parameters governing the magnitude of each possible interaction in the network (Figure 3E). We find that the most important parameters as determined by the sensitivity analysis depends on the choice of sensitivity function. When the sensitivity function for aPAR is used as the model output, we find that aPAR activation by CDC-42 (*q<sup>a</sup>* <sup>3</sup>) is the most important parameter, followed by CDC-42 activation by aPAR (*q<sup>c</sup>* <sup>3</sup>) and pPAR inhibition by aPAR (*k p* <sup>3</sup> ). When the sensitivity function for pPAR is used as the model output, we find that pPAR inhibition by aPAR (*k p* <sup>3</sup> ) is the most important parameter, followed by mutual activation of aPAR and CDC-42 (*q<sup>a</sup>* <sup>3</sup> and *<sup>q</sup><sup>c</sup>* 3) as well as CDC-42 inhibition by pPAR (*k<sup>c</sup>* <sup>3</sup>). Using the sensitivity function for CDC-42 as the model output, we find that CDC-42 activation by aPAR (*q<sup>c</sup>* <sup>3</sup>) is the most important parameter, followed by pPAR inhibition by aPAR (*k p* <sup>3</sup> ), aPAR activation by CDC-42 (*q<sup>a</sup>* 3), and CDC-42 inhibition by pPAR (*k<sup>c</sup>* 3). Combining the most important parameters for each sensitivity function, we find that they correspond to the network interactions in our minimal network (*q<sup>a</sup>* <sup>3</sup>, *<sup>q</sup><sup>c</sup>* <sup>3</sup>, and *k p* <sup>3</sup> interactions, Figure 3B(b3)). Since the minimal network was found independently of the sensitivity analysis, this suggests our minimal network contains only the most critical interactions.

We then evaluated the total-effect index, *STi* . We find that CDC-42 activation by aPAR (*qc* <sup>3</sup> interaction) appears as a high index value with respect to both the aPAR and pPAR sensitivity functions, which we did not find with the first order index (see Appendix A.4 for further details of index value significance and meaning). When considering the aPAR sensitivity function, the mutual activation of aPAR and CDC-42 (*q<sup>a</sup>* <sup>3</sup> and *<sup>q</sup><sup>c</sup>* <sup>3</sup> interactions) have a high total index, indicating that these two mutual activation interactions are very important for the maintenance of aPAR polarity (Figure 3D). Similarly, when considering the pPAR sensitivity function, the two high index values correspond to *q<sup>c</sup>* <sup>3</sup> and *k p* <sup>3</sup> interactions, suggesting that the activation of CDC-42 by aPAR plays an important role in the maintenance of pPAR polarity along with the inhibition of pPAR by aPAR. When considering the sensitivity function for CDC-42, the total-effect index followed the same trend as the first order index, with CDC-42 activation by aPAR (*k<sup>c</sup>* <sup>3</sup> interaction) being the most important parameter. When considering the total effect of a parameter, this suggests that CDC-42 plays a critical role for maintenance of both aPAR and pPAR polarity. We also found that the magnitude of the sensitivity indices for *Si* and *STi* are substantially different for the pPAR and CDC-42 cases, but not for the aPAR case, suggesting that maintenance of pPAR and CDC-42 polarity are more strongly affected by interactions with the other proteins in the network, but maintenance of aPAR polarity operates largely independently of the other proteins. Together, this indicates that aPAR may play a central role for the entire network during maintenance of polarity, and that maintenance of pPAR polarity may depend directly on maintenance of aPAR polarity.

#### *3.2. Critical Network Interactions and Parameters for Generation of Par Protein Polarization*

In the previous section, we determined the critical interactions and parameters in a minimal network for the maintenance of Par protein polarity in the early *C. elegans* embryo. Under wild type conditions, establishment of polarization in the early *C. elegans* embryo relies on actin and myosin based advective flow [17,36]. However, polarization can still occur even in the absence of cortical flow [30], although on a longer time scale. Thus, using a similar approach as in the previous section, we aim to determine if the interactions and parameters required for maintenance of polarity are different from those required to generate polarization in the absence of cortical flow.

#### 3.2.1. Minimal Network for Generation of Par Protein Polarization

To determine the minimal network required to generate a polarity pattern, we simulated networks with individual interactions missing to determine which interactions are predicted to be necessary for self-organization. Simulations were given a local stimulus as an initial conditions, as shown in Figure 2B, first panel, and were assessed at 30 min for the presence of polarization. As shown in Figure 4A, the three interactions represented by the parameters *q<sup>a</sup>* <sup>3</sup>, *<sup>q</sup><sup>c</sup>* <sup>3</sup>, and *k p* <sup>3</sup> play critical roles in generation of polarization as their absence leads to the absence of pattern generation and these are the same three, interactions required for maintenance of polarization (Figure 3B(b3)). However, the minimal network for maintenance of polarization (Figure 3B(b3)) has not been sufficient to generate a polarized pattern (Figure 4B), and the addition of either the *k<sup>a</sup>* <sup>3</sup> or *<sup>k</sup><sup>c</sup>* <sup>3</sup> interaction was also important (Figure 4A(a4,a5)). We have also confirmed that the minimal network for aPAR polarity maintenance (*q<sup>a</sup>* <sup>3</sup> and *<sup>q</sup><sup>c</sup>* <sup>3</sup> interactions, Figure 3D) is not sufficient to generate a polarity pattern even in the presence of CDC-42 inhibition by pPAR (Figure 4C). These results suggest that in addition to the minimal network for maintenance phase, aPAR inhibition by pPAR is critical for generating Par protein polarity, and the mutual inhibition between the anterior and posterior Par proteins directly and via CDC-42 is important in generating the pattern.

**Figure 4.** Minimal network needed for generation of Par protein polarization. (**A**–**C**) Representative simulations for the indicated network. Dotted gray lines indicate interactions excluded from the network. The initial condition provides a small local perturbation, as discussed in Section 2.4. Simulation results are shown at *t* = 30, longer than the time scale for polarity emergence reported experimentally [30]. (**A**) (a1)–(a5) Simulation results of networks omitting a single interaction. Networks boxed in gray are minimal networks capable of generating a polarization pattern. (**B**) and (**C**) (c1)–(c2) Simulation results of networks omitting two or more interactions. (**E**–**F**) Variance-based sensitivity analysis results using the eFAST method. Both the first order and total effect index indicate the importance of CDC-42 in producing the correct spatial pattern.

#### 3.2.2. Critical Parameters for Generation of Par Protein Polarization

We then used sensitivity analysis to determine which parameters are critical for generation of polarization. We used the same sensitivity functions as above (Equation (6)), and assessed the model simulation at *t* ∗ = 30 min. In contrast to our sensitivity analysis results for the maintenance of polarization, we find that CDC-42 activation by aPAR (*q<sup>c</sup>* <sup>3</sup> interaction) is the most important parameter for aPAR, pPAR and CDC-42 with respect to both the first order and total-effect index. Although the interactions represented by the parameters *q<sup>a</sup>* <sup>3</sup>, *<sup>q</sup><sup>c</sup>* <sup>3</sup>, and *k p* <sup>3</sup> tended to show higher importance than the other two interactions, *k<sup>a</sup>* <sup>3</sup> and *<sup>k</sup><sup>c</sup>* <sup>3</sup>, in both first order and total-effect indices, the effect of the *<sup>k</sup><sup>a</sup>* <sup>3</sup> and *<sup>k</sup><sup>c</sup>* 3 interactions are not negligible in the total-effect index as compared to the first order index for aPAR and CDC-42. This indicates that the *k<sup>a</sup>* <sup>3</sup> and *<sup>k</sup><sup>c</sup>* <sup>3</sup> interactions may influence the generation of the polarity pattern, consistent with our minimal network analysis (Figure 4A(a4,a5)).

#### **4. Discussion**

Through minimal network analysis and eFAST sensitivity analysis, we identified different roles for aPAR, pPAR and CDC-42 in maintaining or generating the spatial pattern associated with polarization. To focus on the biochemical interactions, we focus on the dynamics of polarization in the absence of cortical flow, and neglect interactions between CDC-42 and the Par proteins with mechanical proteins including actin and myosin. Results are summarized in Figure 5 and Table 1. CDC-42 is required in both pattern maintenance and generation: CDC-42 reinforces maintenance of aPAR polarity which in turn directs pPAR polarization as indicated by the high sensitivities of the latter to the up-regulation of CDC-42 and down-regulation of pPAR by aPAR. Thus, the entire system of interactions relies on aPAR polarity to maintain polarization in all variables. However, the minimal network capable of maintaining Par protein polarization is insufficient to generate the polarization pattern, and additional inhibition of aPAR or CDC-42 by pPAR is required. This additional interaction acts to balance the mutual inhibition between aPAR and pPAR, allowing a pattern to be generated. In other words, the posterior Par proteins play a more significant role when generating the polarity pattern as compared to maintaining an established pattern. Interactions between CDC-42 and the anterior Par proteins are critical in both situations (Figure 5 and Table 1). This is consistent with observations suggesting the importance of anterior Par proteins and CDC-42 in generating polarization [37]. This suggests that pPAR, with the support of CDC-42, plays a critical role in generating polarization, though, unlike maintenance, we found all interactions between CDC-42, pPAR and aPAR suggested by previous observation, and thus under consideration, can have some influence on polarization generation.

**Figure 5.** Critical interactions during maintenance and generation of Par protein polarization. A bigger circle indicates a key role for that protein in the indicated phase. In the maintenance phase, aPAR plays a key role in maintaining spatial polarity via interactions with CDC-42. aPAR enforces pPAR polarity through mutual inhibition. In the patterning phase, pPAR is playing the key role by inhibiting aPAR through either direct or CDC-42-mediated inhibition. CDC-42 is indispensable for both the maintenance and emergence of polarization, but does not play a critical role in the dynamics of either phase.

Remarkably, the minimal network analysis and eFAST both selected the same critical interactions and parameters in polarization maintenance and generation. This independent determination suggests the minimal networks found here represents the key interactions involved in patterning associated with polarization.

**Table 1.** Summary of minimal network and sensitivity analyses results. MNA: Minimal Network Analysis, *Si*: First order index, *STi* : Total effect index. For MNA, a small checkmark indicates additional interactions required for pattern generation. For sensitivity measures, a/p/c denote the aPAR/pPAR/CDC-42 sensitivity functions, respectively, and a capital letter indicates the parameter with the highest index for that function.


In this study, we have made a number of simplifying assumptions to investigate the core interactions between the Rho protein CDC-42 and members of the Par protein family. In future investigations, we will explicitly consider separate members of the Par protein family, each of which have distinct dynamics. For instance, the anterior Par protein PAR-3 appears to compete with CDC-42 to bind to a complex of PAR-6 and aPKC [18], and in the posterior, PAR-2 must be recruited to the membrane before PAR-1 can bind [10]. We will add other members of the Rho protein family, including Rac and Rho, which are known to act in a network with CDC-42 [38]. In this investigation, we have explored protein interactions within the pattern maintenance phase of early embryo development, when advective flow has largely ceased [17,22], and in addition considered the differences in protein interactions required for pattern generation in the absence of cortical flow, as observed in select experimental systems. Future extensions of the model will include the spatial dynamics of advective flow, which may act with the biochemical interaction network to establish polarity. In this study, we have reduced the geometry of the embryo to a 1D domain for computational efficiency to facilitate the exploration of a five dimensional parameter space. However, other studies have studied polarization dynamics of simplified two variable models on fully 3D domains [39] or explicitly included the geometry of the embryo [40] to investigate the spatial orientation of the polarity axis. Noting the computational demands, we leave minimal network and eFAST analysis of higher dimensional parameter spaces in more complex geometries and 3D domains for future work.

In this investigation, we have developed and analysed a mathematical model of the biochemical network that integrates components of two protein families, Par proteins and Rho GTPases, that are known polarity determinants. The minimal network and sensitivity analyses have identified critical interactions in this network, providing testable hypotheses for future experimental work. By resolving the critical role of CDC-42 and highlighting the most important PAR interactions in this network, we have extended our understanding of the signalling network responsible for polarization, and laid the foundation for further investigations into patterning associated with polarization.

**Author Contributions:** Conceptualization, S.S.-L., E.A.G., A.T.D.; methodology, S.S.-L., E.A.G., A.T.D.; investigation, S.S.-L., E.A.G., A.T.D.; writing—original draft preparation, S.S.-L., E.A.G., A.T.D.; writing—review and editing, S.S.-L., E.A.G., A.T.D. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by a Grant-in-Aid for Scientific Research from the Ministry of Education, Culture, Sports, Science and Technology, Japan, to S.S.L. (JP19H01805 and JP17KK0094) and funding from the National Science Foundation (USA) to A.T.D. (DMS-1554896).

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **Appendix A**

#### *Appendix A.1. Model Reduction*

The model reduction is summarized schematically in Figure 2A; here we provide the associated simplifications of the governing equations, starting from Equation (3), that is

$$\begin{cases} \frac{\partial [A\_m]}{\partial t} = D\_m^A \nabla\_\Gamma^2 [A\_m] + \left\{ \gamma\_d + F\_{\text{off}}^A([\mathbb{C}\_m]) \right\} \left( \frac{A\_{\text{tot}}}{|\Omega|} - \frac{1}{|\Omega|} \int\_\Gamma [A\_m] d\mathbf{x} \right) - \left\{ a\_d + F\_{\text{off}}^A([P\_\mathcal{U}]) \right\} [A\_m] & \text{ on } \mathbf{x} \in \Gamma, \\\frac{\partial [P\_\mathcal{W}]}{\partial t} = D\_m^P \nabla\_\Gamma^2 [P\_\mathcal{W}] + \gamma\_p \left( \frac{P\_{\text{tot}}}{|\Omega|} - \frac{1}{|\Omega|} \int\_\Gamma [P\_m] d\mathbf{x} \right) - \left\{ a\_p + F\_{\text{off}}^P([A\_m]) \right\} [P\_\mathcal{W}] & \text{ on } \mathbf{x} \in \Gamma, \\\frac{\partial [\mathcal{C}\_m]}{\partial t} = D\_m^C \nabla\_\Gamma^2 [\mathcal{C}\_m] + \gamma\_\mathcal{c} \left( \frac{\mathcal{C}\_{\text{tot}}}{|\Omega|} - \frac{1}{|\Omega|} \int\_\Gamma [\mathcal{C}\_m] d\mathbf{x} \right) - \left\{ a\_c + F\_{\text{off}}^C([P\_\mathcal{W}], [A\_m]) \right\} [\mathcal{C}\_m] & \text{ on } \mathbf{x} \in \Gamma. \end{cases} \tag{A1}$$

Axisymmetry of the cell geometry is assumed about the major axis of symmetry of the cell in Figure 2A, and the objective is to further simplify model (3) to the one dimensional model of (4). To proceed, let *R*(*s*) denote the radial distance of the membrane from the major axis of symmetry, with *s* ∈ [0, *L*] the arclength along the perimeter of the depicted, projected membrane in Figure 2A and *L* the total perimeter length. In addition, we take *s* = 0 to be the leftmost point on the membrane, and impose periodic boundary conditions at *s* = 0, *L*. Note that the metric tensor components for this surface with respect to coordinates *X*<sup>1</sup> = *s* and *X*<sup>2</sup> = *φ*, the angle of rotation about the major axis of symmetry, are given by *gij* = diag(1, *R*2(*s*))*ij*. Hence the Laplacian operation for an axisymmetric membrane concentration, denoted *um*(*s*) in general, is given by

$$\nabla\_{\Gamma}^{2}u\_{\mathrm{ff}} = \frac{1}{\mathcal{S}^{1/2}}\frac{\partial}{\partial X\_{i}}\left(\mathcal{g}^{1/2}\mathcal{g}^{j}\frac{\partial u\_{\mathrm{m}}}{\partial X\_{j}}\right) = \frac{1}{R(s)}\frac{\partial}{\partial s}\left(R(s)\frac{\partial u\_{\mathrm{m}}}{\partial s}\right) = \frac{\partial^{2}u\_{\mathrm{m}}}{\partial s^{2}} + \frac{R\_{s}}{R}\frac{\partial u\_{\mathrm{m}}}{\partial s},\tag{A2}$$

where the subscript s denotes the ordinary derivative with respect to arclength, summation convention is used, *g* = det(*gij*), *gij* = (*gij*)−<sup>1</sup> and noting that *φ*-derivatives generate no contribution by axisymmetry.

We also note that the mechanism of patterning for the above equations is widely recognized to be wave pinning, due to the multiple steady states associated with the kinetics. The pattern organization that emerges corresponds to the emergence of a wave transitioning between different steady states, with homogeneous concentrations elsewhere [41]. In particular, the location where the transition wave halts is approximated by an estimate of where the speed for wave-like solutions of the governing equations drops to zero, so that the final pattern is a standing wave, rather than, for instance, a Turing pattern [41].

Furthermore, for the presented model, these transitions between steady states are sharp for the parameter values of Table A1, as also seen *a posteriori* in Figure 2B for example. This is a consequence of the small size of the *non-dimensional* diffusion coefficients relative to the *non-dimensional* inverse timescales in Table A1, with <sup>2</sup> ∼ *<sup>D</sup><sup>A</sup> <sup>m</sup>*/*γ<sup>a</sup>* ∼ <sup>1</sup>×10−3, so that the lengthscale of these transitions corresponds to an order unity coefficient multiplying 1.

Thus either side of a transition is a homogeneous solution, whereby *um* is at a constant steady state for each concentration in the model and the advective term [*Rs*/*R*]*∂um*/*∂s* is not important. Within a transition, there are steep gradients generating a dominant balance between the Laplacian and the kinetics. In particular for any given concentration, after non-dimensionalization and with the non-dimensional kinetics generically denoted by *K*∗, we have equations of the form

$$\eta \frac{\partial \mu\_{\ast \text{m}}}{\partial t\_{\ast}} = \epsilon^2 \nabla\_{\ast \Gamma}^2 \mu\_{\ast \text{m}} + K\_{\ast} = \epsilon^2 \left( \frac{\partial^2 \mu\_{\ast \text{m}}}{\partial s\_{\ast}^2} + \frac{R\_{\ast s\_{\ast}}}{R\_{\ast}} \frac{\partial \mu\_{\ast \text{m}}}{\partial s\_{\ast}} \right) + K\_{\ast \text{s}} \tag{A3}$$

where asterisks indicate non-dimensionalized variables and *η* is a non-dimensional constant that emerges from the rescaling. With the further change of variable *S* = (*s*<sup>∗</sup> − *s*0)/, where *s*<sup>0</sup> is the location of the transition, to make the dominant balance explicit by matching the transport and kinetic terms and to generate order unity derivatives, we have that the inner region equations for the transition layer in a matched asymptotic expansion approximation are given by

$$
\eta \frac{\partial u\_{\ast \mathfrak{m}}}{\partial t\_{\ast}} = \frac{\partial^2 u\_{\ast \mathfrak{m}}}{\partial S^2} + \epsilon \left( \left. \frac{R\_{\ast \mathfrak{s}\_{\ast}}}{R\_{\ast}} \right|\_{s\_0} + O(\epsilon S) \right) \frac{\partial u\_{\ast \mathfrak{m}}}{\partial S} + K\_{\ast} = \frac{\partial^2 u\_{\ast \mathfrak{m}}}{\partial S^2} + K\_{\ast} + O(\epsilon), \tag{A4}
$$

with the last step valid provided *R*∗*s*∗/*R*<sup>∗</sup> ∼ *O*(1/). Under such circumstances the advection term does not contribute to the structure of the equations either in the outer regions where the solution is at steady state, or the inner region. In turn, this justifies approximating the the Laplacian with the second derivative ∇<sup>2</sup> <sup>Γ</sup>*um* ≈ *<sup>∂</sup>*2*um*/*∂s*2. Furthermore, with this scaling, where all terms are order unity, it can be seen that the advective term will only generate a subleading contribution to the estimates of where the transitions become fixed as standing waves [41], so that the temporal dynamics of interest, that fix the location of the pattern, are also well approximated in the absence of the advective term away from regions where *R*∗*s*∗/*R*<sup>∗</sup> ∼ *O*(1/).

However, sufficiently near the poles at *s* = 0, *L*/2, 1/*R* becomes arbitrarily large. With *X* denoting distance along the horizontal axis of symmetry in Figure 2A, we have *s*<sup>2</sup> *<sup>R</sup>* = <sup>1</sup> + 1/*R*<sup>2</sup> *<sup>X</sup>* and, noting *RX* blows upon approaching the poles, we have *R*(*s*) ≈ *s* and, in turn *R*<sup>∗</sup> ≈ *s*∗, for *s*<sup>∗</sup> . Thus *R*∗*s*∗/*R*<sup>∗</sup> ≈ 1/*s*<sup>∗</sup> ∼ *O*(1/) once *s*<sup>∗</sup> ∼ *O*() near the left hand pole and analogously near the right hand pole. Thus, approximating the Laplacian by second derivatives only breaks down within a distance of of the poles in the non-dimensional model, and only if there is a transition in this region. In particular, with the appropriate scaling of *S* = *s*∗/ for a potential transition within of the left hand pole, we have

$$
\eta \frac{\partial u\_{\ast m}}{\partial t\_{\ast}} = \frac{\partial^2 u\_{\ast m}}{\partial S^2} + \frac{1}{S} \frac{\partial u\_{\ast m}}{\partial S} + K\_{\ast \prime} \tag{A5}
$$

and the equations take the form of the well-studied, radially symmetric, system near the origin of the 2D plane. The impact of the advective term here is to slow the transition wave, though its wavespeed asymptotes to that associated with Equation (A4) once the transition has propagated into a region with *s*<sup>∗</sup> (e.g., Eqn. 11.20 *et. seq.* in Murray's textbook (1993)), with analogous conclusions for the right hand pole at *s* = *Ls*<sup>∗</sup> = *L*/2.

In the presented model we only consider the second derivatives. Then with patterns emerging from one of the poles (for example as in Figure 2B, with an emerging pattern at *x* ≡ *s* = *Ls*<sup>∗</sup> = *L*/2), the presented model will underestimate the timescale for the emergence of the wave, but otherwise will capture its dynamics once the transition is beyond a distance of from *x* = *L*/2. Furthermore, unless a steady state transition occurs -close to the left or right poles of the cell, its location and spatial variation will be well approximated even in the absence of the advective term in the model equations, as motivated above. Thus, and in summary, in the simplified model of (4) we take ∇<sup>2</sup> <sup>Γ</sup>*um* = *∂*2*um*/*∂s*<sup>2</sup> for all concentrations, assuming that the membrane changes sufficiently slowly in shape to ensure *R*∗*s*∗/*R*<sup>∗</sup> = *LRs*/*R* ∼ *O*(1/) holds away from the poles. Then apart from the temporal dynamics as a transition passes through, or emerges from, a pole the simplified model yields an accurate picture, especially for the final steady-state pattern, providing a transition does not occur at the poles (and a posteriori the results presented do not indicate steady state transitions at the poles).

The approximation in particular greatly simplifies the model in that it sidesteps dealing with coordinate singularities at the poles. While these are possible to accommodate, this entails additional numerical complexity, such as Taylor expanding and asymptotically approximating near the poles (Morton and Mayers [42], Woolley et al. [43]) or, when such approximations are not feasible, borrowing from basic differential geometry and working with an atlas of charts and the associated multiple coordinate systems [44], which in this model would result in the loss of the simplifications

from axisymmetry. Such complexities are not warranted for the current study given the need for extensive computation with the global sensitivity analyses and the minor impact of the approximation ∇2 <sup>Γ</sup>*um* ≈ *<sup>∂</sup>*2*um*/*∂s*2, where *<sup>x</sup>* ≡ *<sup>s</sup>* in the main text.

We further approximate the geometry by a slender cylinder radius *H*, length L ≈ *L*/2 − 2*H*, with essentially neglected asymptotically short caps of length *δ* with *δ H L*. Then, on use of axisymmetry, the surface integrals within the governing Equation (3) simplify via

$$\begin{split} \int\_{\Gamma} u\_{\mathrm{m}} \mathrm{d}x &= \,^{2}\pi \int\_{0}^{L/2} u\_{\mathrm{m}} \mathrm{R} \mathrm{d}x = 2\pi H \int\_{0}^{L/2} u\_{\mathrm{m}} \mathrm{d}x + 2\pi \int\_{0}^{P} u\_{\mathrm{m}} (\mathrm{R} - H) \mathrm{d}x + 2\pi \int\_{L/2 - P}^{L/2} u\_{\mathrm{m}} (\mathrm{R} - H) \mathrm{d}x \\ &= \,^{2}\pi H \int\_{0}^{L/2} u\_{\mathrm{m}} \mathrm{d}x \left( 1 + O\left(\frac{H}{L}\right) \right) = \int\_{0}^{L} \pi H u\_{\mathrm{m}} \mathrm{d}x \left( 1 + O\left(\frac{H}{L}\right) \right), \end{split} \tag{A6}$$

where *x* is dimensional arclength, *P* ≈ *H* is the contribution to the cell perimeter arclength from the left hand pole to the first point where *R* = *H*. Noting *R* ≈ *x* for the left hand cap and similarly for the right hand cap the integrals can be approximated as above to accuracy *H*/*L* given the separation of the geometrical scales. These corrections are also dropped in the final dimensional model given *<sup>H</sup> <sup>L</sup>*. Finally it is easier to work with [*A*˜*m*] = *<sup>π</sup>H*[*Am*], which is the density per unit length away from the caps, and we analogously define [*P*˜ *<sup>m</sup>*], [*C*˜*m*]. Then

$$\begin{split} \frac{\partial[\bar{A}\_{m}]}{\partial t} &= D\_{m}^{A} \frac{\partial^{2}}{\partial \mathbf{x}^{2}} [\bar{A}\_{m}] + \left\{ \bar{\gamma}\_{\text{f}} + F\_{\text{en}}^{A} ([\bar{\mathsf{C}}\_{m}]) \right\} \left( \frac{A\_{\text{tot}}}{L} - \frac{1}{L} \int\_{0}^{L} [\bar{A}\_{m}] dx \right) - \left\{ a\_{\text{t}} + F\_{\text{off}}^{A} ([P\_{\text{m}}]) \right\} [\bar{A}\_{m}], \\ \frac{\partial[\bar{P}\_{m}]}{\partial t} &= D\_{m}^{\overline{\mathsf{D}}} \frac{\partial^{2}}{\partial \mathbf{x}^{2}} [\bar{P}\_{m}] + \bar{\gamma}\_{\text{P}} \left( \frac{P\_{\text{tot}}}{L} - \frac{1}{L} \int\_{0}^{L} [\bar{P}\_{\text{m}}] dx \right) - \left\{ a\_{\text{p}} + F\_{\text{off}}^{P} ([\bar{A}\_{m}]) \right\} [\bar{P}\_{\text{m}}], \\ \frac{\partial[\bar{\mathcal{C}}\_{m}]}{\partial t} &= D\_{m}^{\overline{\mathsf{C}}} \frac{\partial^{2}}{\partial \mathbf{x}^{2}} [\bar{\mathsf{C}}\_{m}] + \bar{\gamma}\_{\text{c}} \left( \frac{\mathsf{C}\_{\text{tot}}}{L} - \frac{1}{L} \int\_{0}^{L} [\mathcal{C}\_{m}] dx \right) - \left\{ a\_{\text{c}} + F\_{\text{off}}^{C} ([\bar{P}\_{\text{m}}], [\bar{A}\_{m}]) \right\} [\bar{\mathsf{C}}\_{m}], \end{split} \tag{A7}$$

where *<sup>γ</sup>*˜ *<sup>a</sup>* = *<sup>γ</sup>aLπH*/|Ω|, with analogous definitions for *<sup>γ</sup>*˜ *<sup>p</sup>*, *<sup>γ</sup>*˜ *<sup>c</sup>* and also for *<sup>q</sup>*˜*<sup>a</sup>* <sup>3</sup> = *<sup>q</sup><sup>a</sup>* <sup>3</sup>*πHL*/|Ω|. Dropping tildes, and thus redefining [*Am*], [*Pm*], [*Cm*], *γa*, *γp*, *γc*, *q<sup>a</sup>* <sup>3</sup>, gives the final dimensional model, Equation (4) of Section 2.2.

#### *Appendix A.2. Additional Simulations and Representative Parameter Set*

The simulation results for Model (1) with fast cytosol diffusions are shown in Figure A1. The representative parameter sets which we used in simulations are given in Table A1.

**Figure A1.** Simulations results for Model (1) on 1D. (**A**) The simulation results of Model (1) with biologically feasible cytosol diffusion coefficient. (**B**) The simulation results for the case where the cytosol diffusions are ten times larger than the case of (**A**). Both cases (**A**) and (**B**) show almost homogeneous concentrations for cytosol proteins and with a sufficiently large diffusion coefficient, we can assume that the asymptotic dynamics of cytosol proteins is homogeneous steady state. The simulation results are reported at 30 min and the detailed parameters are given in Table A1.

**Table A1.** Representative parameter set. The parameters *L*, *τ*, M are respectively a lengthscale, timescale and a scale of protein number used to non-dimensionalize the model. The dimensional parameters are defined in terms of the protein number scale, so that the non-dimensional model is independent of M. In practice, this means the parameters associated with the non-linear kinetics are such that the non-linear reactions are significant and also in balance with other reactions at typical cellular numbers of proteins. It should be noted that the parameter set below is inherently non-identifiable as, for example, the model depends on the three parameters *q<sup>a</sup>* <sup>1</sup>, *<sup>q</sup><sup>a</sup>* <sup>2</sup>, *<sup>q</sup><sup>a</sup>* <sup>3</sup> only via the two degrees of freedom *q<sup>a</sup>* 1/*q<sup>a</sup>* <sup>3</sup>, *<sup>q</sup><sup>a</sup>* 2/*q<sup>a</sup>* <sup>3</sup>. For presentational simplicity we have not eliminated such superfluous degrees of freedom from the parameter set, though this would be necessary if parameter inference from experimental data was pursued in future work.


†1 [30], †2 [34], †3 [45].

*Appendix A.3. Variance-Based Sensitivity Analysis by Using Extended Fourier Amplitude Sensitivity Test (eFAST) Method*

For a large enough sample size, the variance-based sensitivity analysis provides a quantitative measure for how much factor *A* is more important than factor *B*. It generally requires extensive computation and the FAST method has been demonstrated as a method to reduce the cost effectively by exploring the multidimensional space of the input factors via a suitably defined search-curve [31,32]. eFAST is a generalization the FAST method that provides the total effect index, defined below.

Let us define *<sup>y</sup>* as the expected model output and *pi*(*<sup>i</sup>* <sup>=</sup> <sup>1</sup> ··· *<sup>N</sup>*¯ ) are the input factors (model parameters). The effect (contribution) of each parameter *pi* is considered as the variance of the expected model output, namely, *Vi* = *V*(*E*(*y*|*pi*)). The effect of the interaction between two orthogonal inputs *pi* and *pj* on the output *y* is defined in terms of conditional variances as *Vij* = *V*(*E*(*y*|*pi*, *pj*)) − *Vi* − *Vj* (Sections 5.9 and 5.10, [31]). Similarly, the effect of the interaction among three distinct inputs *pi*, *pj* and *pk* on the output *y* is given as *Vijk* = *V*(*E*(*y*|*pi*, *pj*, *pk*)) − *Vij* − *Vik* − *Vjk* − *Vi* − *Vj* − *Vk* and so on. Then the total output variance *Vtotal* for a model with *N*¯ input factors is given by

$$V\_{total} = \sum\_{i} V\_i + \sum\_{i} \sum\_{j>i} V\_{ij} + \sum\_{i} \sum\_{j>i} \sum\_{k>j} V\_{ijk} + \dots + V\_{1\dots N}.$$

The first order sensitivity index is defined by

$$S\_i = \frac{V\_i}{V\_{total}}\_{\prime}$$

and the total-effect index is defined by

$$S\_{T\_i} = \frac{V\_{total} - V(E(y|p\_{-i}))}{V\_{total}} = 1 - \frac{V\_{-i}}{V\_{total}}e$$

where −*i* stands for *all but i* and *V*−*<sup>i</sup>* is given to

$$V\_{-i} = V\_{total} - V\_i - \sum\_{j} V\_{ij} - \dots - V\_{1\dots \hat{N}} \cdot 1$$

The sensitivity analysis using eFAST has been carried out according to the following steps [32,35].

#### **STEP 1: Sampling for search-curve**

Each parameter value was sampled according to

$$p\_i(s) = p\_i^{\min} + (p\_i^{\max} - p\_i^{\min}) \left[ \frac{1}{2} + \frac{1}{\pi} \sin^{-1} (\sin(\omega\_i s + \psi\_i)) \right],$$

where *s* ∈ (−*π*, *π*) and *ψ<sup>i</sup>* is a random phase-shift chosen in [0, 2*π*) which is used for making different curves (i.e., resampling) and we carried out two different curves by resampling. {*ωi*} is a set of different angular frequencies associated with each input factor. The sample size (*Ns*) is given by *Ns* = (2*Mω*max + 1)*Nr* where *M* is the interference factor (usually 4 or higher), *ω*max is highest frequency and *Nr* is the number of resamplings. The detailed values which we chosen are shown in Table A2.

**Table A2.** The values used for sensitivity analysis.


#### **STEP 2: Calculating output of sensitivity function**

With respect to the parameter samples, we numerically solved the main model (4) and calculated each sensitivity function, **y**(*s*) = *FS*(**p**(*s*)), given by the Equation (6), where **p**(*s*)=(*p*1(*s*), ··· , *pN*¯ (*s*)) and **y** = (*y*1(*s*), ··· , *yN*¯ (*s*)), with *s* partitioned via the nodes *sk* = *π*(2*k* − *Ns* − 1)/*Ns*, *k* ∈ {1, ··· , *Ns*}.

#### **STEP 3: Calculating the power spectrum in Fourier series**

We expand *FS*(*s*) in a Fourier series such that

$$F\_S(s) = \sum\_{j = -\infty}^{+\infty} \{A\_j \cos(js) + B\_j \sin(js)\},$$

where the Fourier coefficients *Aj* and *Bj* are defined as

$$\begin{aligned} A\_{\circ} &= \frac{1}{2\pi} \int\_{-\pi}^{\pi} F\_{\mathsf{S}}(s) \cos(js) ds = \frac{1}{2\pi} \sum\_{k=1}^{Ns} F\_{\mathsf{S}}(s\_k) \cos(js\_k) \triangle s = \frac{1}{N\_s} \sum\_{k=1}^{Ns} F\_{\mathsf{S}}(s\_k) \cos(js\_k), \\ B\_{\circ} &= \frac{1}{2\pi} \int\_{-\pi}^{\pi} F\_{\mathsf{S}}(s) \sin(js) ds = \frac{1}{2\pi} \sum\_{k=1}^{Ns} F\_{\mathsf{S}}(s\_k) \sin(js\_k) \triangle s = \frac{1}{N\_s} \sum\_{k=1}^{Ns} F\_{\mathsf{S}}(s\_k) \sin(js\_k), \end{aligned}$$

and *s* = 2*π*/*Ns*. Then the power spectrum is calculated by

$$
\Lambda^2(j) = A\_j^2 + B\_j^2.
$$

#### **STEP 4: Calculating** *Si* **and** *STi*

The total variance (*Vtotal*), the variance of *i* factor (*Vi*), and the variance of *all but i* (*V*−*i*) are calculated by

$$\begin{aligned} V\_{\text{total}} &\approx \mathcal{V}\_{\text{total}} = 2 \sum\_{\omega=1}^{M\omega\_{\text{max}}} \Lambda^2(\omega), \\ V\_i &\approx \mathcal{V}\_i = 2 \sum\_{p=1}^M \Lambda^2(p\omega\_i), \\ V\_{-i} &\approx \mathcal{V}\_{-i} = 2 \sum\_{p=1}^{\frac{\omega\_i}{2}} \Lambda^2(p\omega), \end{aligned}$$

as in [35]. We average the variances above over the sets of resampling and calculate the indices

$$S\_i = \frac{V\_i}{V\_{total}} \qquad \text{and} \qquad S\_{T\_i} = 1 - \frac{V\_{-i}}{V\_{total}}.$$

Note that we need to give the maximal frequency value for *ω<sup>i</sup>* of *i th* factor in calculating *<sup>V</sup>*−*i*, so that total-effect variance is calculated separately for each factor *i*.

#### *Appendix A.4. First Order Index and Total Index for the Dummy Parameter*

Since the eFAST method artifactually produces small but non-zero sensitivity indices, we have calculated the first order index of the dummy parameter (*S*dummy) and the total index of the dummy parameter (*ST*dummy ) in order to confirm that the index values which we obtained in Figures 3 and 4 give reliable data to see the significance of model parameters [46]. If the model parameters with a total index less than or equal to that of the dummy parameter, the model parameters should be considered not significantly different from zero.

For the calculations of the first order index for the dummy parameter, we put the frequency of the dummy parameter to *ω*dummy = 91 with the same frequency set given in Table A2. In the calculation of the total index for the dummy parameter, we used additional frequency *ω* = 13 together with the same frequency set given in Table A2. The sensitivity results for the dummy parameter are shown in Table A3.


**Table A3.** First order index and total index for the dummy parameter.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
