Cytoskeleton polarity is essential in determining orientational order in basal bodies of multi-ciliated cells.

In multi-ciliated cells, directed and synchronous ciliary beating in the apical membrane occurs through appropriate configuration of basal bodies (BBs, roots of cilia). Although it has been experimentally shown that the position and orientation of BBs are coordinated by apical cytoskeletons (CSKs), such as microtubules (MTs), and planar cell polarity (PCP), the underlying mechanism for achieving the patterning of BBs is not yet understood. In this study, we propose that polarity in bundles of apical MTs play a crucial role in the patterning of BBs. First, the necessity of the polarity was discussed by theoretical consideration on the symmetry of the system. The existence of the polarity was investigated by measuring relative angles between the MTs and BBs using published experimental data. Next, a mathematical model for BB patterning was derived by combining the polarity and self-organizational ability of CSKs. In the model, BBs were treated as finite-size particles in the medium of CSKs and excluded volume effects between BBs and CSKs were taken into account. The model reproduces the various experimental observations, including normal and drug-treated phenotypes. Our model with polarity provides a coherent and testable mechanism for apical BB pattern formation. We have also discussed the implication of our study on cell chirality.

Introduction

Synchronous and coordinated ciliary beating underlies various biological functions in organisms ranging from ciliates to mammals [1]. In the mammalian trachea, synchronous ciliary beating of multi-ciliated cells produces mucociliary transport, that is the flow of mucus to the oral side of the epithelial cell surface, ejecting viruses and dust to the outside of the body [2]. For efficient mucociliary transport, hundreds of cilia covering the apical membrane of a cell beat synchronously in the same direction. This organized directionality in ciliary beating is achieved by basal bodies (BBs), which are roots of the cilia. These BBs point in the same orientation (Fig 1A), which is defined by the relative position of the basal foot (BF), a major appendage associated with BBs. In addition, previous studies have revealed that BBs are regularly aligned in the apical membrane of the mature epithelial cells in mouse trachea [3]. The positional and orientational order of BBs is crucial for proper ciliary beating, as reported in many other systems where disturbance in ciliary organization resulted in dysfunction of organs, such as inner ear hearing loss [4], infertility [5], hydrocephalus [6], and situs inversus [7]. Therefore, elucidation of the mechanism of ciliary organization on the apical membrane is of paramount biological and clinical importance. However, it is largely unknown how BBs are coordinated through the maturation of the tracheal multi-ciliated cells.

Fig 1

Schematics of the positional and orientational orders of BBs in the trachea.

(A) The mammalian trachea sends mucus to the mouth on the surface of tracheal epithelial cells (red arrow). Tracheal epithelia is composed of multi-ciliated (green) and mucous cells (red). The mucus flow is produced by beating of multiple cilia in the apical membrane of the multi-ciliated cells. The orientation and positioning of cilia are determined by BF (red corn) and BBs (green rod), respectively. BF is an MT-organizing center [8,9]. (B) Typical pattern of BB distribution. During maturation of multi-ciliated cells, the pattern of BBs is changed from cluster to alignment patterns [3]. (C) Polar MTs can orient BBs to the definite direction. (i) Binding of the BF to the neighboring MT bundles is not sufficient to realize unidirectional orientation of BBs [10]. (ii) BBs can orient to same direction if polarity exists in the MT bundles (denoted as ) and BBs can recognize it. The polarity is also shown by blue arrows in panel A. (iii) TEM images of BBs and detected MTs (adapted from [3]). BBs contact MTs at a biased angle. (D) Histograms of the angle of BFs measured using TEM (From (I) Movie S8 in [11], (II) Fig 6 in [11] (III) Video 6 in [3], and (IV) Video 7 in [3]).

Long-term live cell imaging of cultured tracheal cells revealed that maturation of BBs in the apical membrane proceeds though several stages, characterized by a typical pattern in BB distribution (Fig 1B) [3]. Clusters of BBs are observed in the early stages after differentiation to multi-ciliated cells. The clusters are then scattered, followed by formation of short and branched arrays of several BBs termed as ‘partial alignment’ [3]. Through connection and reconnection of these arrays, BBs finally become regularly aligned over the apical membrane. The number of BBs remains nearly unchanged during the maturation period. BBs become increasingly oriented in the same direction indicated by the relative position of the BF, concurrent with the alignment of BBs. Quantification of positional and orientational orders of BBs showed that these two ordering processes in position and orientation were highly correlated [3], suggesting that both are coordinated by common factors. Cytoskeletons (CSK) in the apical membrane were shown to be involved in the regulation of BB distribution. During maturation of cells, microtubules (MTs) and 10 nm-filaments, which are sparsely distributed in the early stage, grow to a dense mesh network and surround the BBs [12]. Most MT bundles run in parallel between the aligned BB array, which might contribute to the positioning formation of the BB merely by invading and separating the array of BBs [11]. Further, since BF directly binds to MTs as confirmed by transmission electron microscopy (TEM) for which γ-tubulin in BF might mediate the association [9], it was suggested that the orientation of the BBs is also regulated by MTs [11]. Remarkably, inhibition of MT polymerization by nocodazole treatment results in disruption of BB alignment pattern and formation of BB cluster pattern, with disruption in the orientation of BB as well. Removal of nocodazole by washing results in recovery of BB positional order [3]. In summary, there is sufficient experimental evidence to show that CSKs are responsible for the determination of position as well as orientation of BBs. However, the underlying mechanism explaining these observations remains elusive. In particular, the attachment of BF to MTs might be important for the determination of BB orientation [9], but its mechanistic details have not been identified as yet.

In the present study, we have generated a mathematical model to explain the ordering of position and orientation in BBs. Our aim here is not to propose a detailed model that can explain all aspects of the experimental findings, but rather to provide a possible model to explore the role of MTs for the observed BB patterning. Firstly, the symmetry of the system has been taken into consideration to show that determination of orientation in BBs requires polarity in MT bundles. The presence of polarity is supported by a quantitative relationship between BB orientation and CSK bundles, determined from available experimental data [3,11]. This was followed by derivation of a mathematical model based on the characteristics of CSKs and their interaction with BBs. ‘Active matter’ formalism was employed by which the generic form of system equations is obtained by taking into account the symmetry of the system [13,14]. This kind of formulation was conducted for studying pattern formation on apical membranes for a system without BBs [15]. Our early study on tracheal cells also employed this scheme and produced positional order of BBs [3], but did not consider the orientation of BBs and the CSK-BB interaction. Using the derived model, we have investigated the conditions for BBs to align and point in the same direction, which further confirms that MT polarity is crucial for the establishment of BB pattern. The information transmission from planar cell polarity (PCP) to the BB orientation has also been discussed. In the Discussion, we have summarized our findings and challenges that need to be addressed in future. We have also considered the implication of our findings in determining a possible mechanism for cellular chirality in determining the relationship between PCP and CSK polarity.

Results

Polarity in cytoskeleton guides orientation of BBs

During the maturation period of multi-ciliated cells, BBs establish regularity in their position and orientation. The patterning of BB arrays follows the formation of the CSK matrix, particularly the MT bundles forming on the cytosolic apical membrane [3,11]. The BBs orient in the same direction at the same time [3]. Since the BF, which is an appendage of BBs, has a component binding to MTs as confirmed by TEM images and mutant strains in Odf2 [8,9,11], it was suggested that the connection of BF with neighboring MT bundles is a determinant of BB orientation [10]. However, this binding between BF and MTs allows BBs to associate with either of the two neighboring MT bundles and hence they cannot distinguish between the two opposite directions (Fig 1C (i)). Therefore, how BBs can determine their definite direction remains a mystery.

A probable hypothesis to fill this knowledge gap is the presence of polarity in the MT bundles. It may be assumed that there is a component in CSK that provides the direction which the BBs are able to interpret. The presence of polarity therefore enables BBs to distinguish the direction by breaking the symmetry in the relationship between BBs and MT bundles (Fig 1C (ii)). To assess this hypothesis, we focused on the contact angles between BBs and MT bundles, since the angles are apparently neither perpendicular nor uniformly distributed but are specifically biased. We measured the angles between MTs and the orientation of the associated BBs from published TEM images (Fig 1C (iii)) [3,11]. The contact angles were biased to the positive side as shown in Fig 1D and Table 1 (n = 4). This observation indicates a polar nature in MTs and chirality in BF, which enables the connection of BBs to MTs within a specific range of angles.

Table 1

Measured angles of BFs and its numbers in Fig 1D.

TEM images	mean	S.D.	numbers	Reference
I	47.2°	±10.5°	10	Movie S8 in [11]
II	41.2°	±8.6°	9	Fig 6 in [11]
III	26.7°	±20.6°	23	Video 6 in [3]
IV	24.3°	±19.0°	36	Video 7 in [3]

Mean and standard deviation of contact angle of BBs to MTs. Each sample is adopted from TEM images shown in Reference. Mean and standard deviation are defined according to directional statistics [16].

Mathematical model for the interaction between CSK and BB

A mathematical model was generated based on the polarity of MTs, demonstrating that polarity is a necessary and sufficient condition for the reproduction of the observed regular positional and orientational orders of BBs on the apical membrane. The following processes were taken into account while generating the model (Fig 2);

Fig 2

Model for positional and orientational orders of BBs.

(A) Illustration of the processes considered in the model. (i) CSK filaments are condensed and bundled in parallel and their polarity orient in the same direction (blue arrow, ). (ii, iii) One BB repels the CSK bundle and other BBs. (iv) Polymerization and depolymerization of CSKs. (v) The orientation of BB, θ, is determined by CSK polarity. (B) Typical time series of BB pattern development (α = 2.0). Green indicates CSK concentration, and blue arrow indicates the polarity of CSK. Gray disk is a BB and the line segment in the disk indicates the orientation of BB. Patten of BBs reaches “alignment”. (C) Time series when the coupling between polarity and CSK concentration is weak (α = 0.2). Eventual pattern is “partial alignment”. Bars, 0.5 μm.

CSKs form bundles on the apical membrane that sustain the intrinsic polarity.

CSKs and BBs exclude each other.

BBs repel each other.

Polymerization and depolymerization of CSKs.

BBs orient in the direction determined by CSK polarity.

The processes (i)-(iii) are modeled by a free energy function F, a generic form of which is as follows [13-15,17].

Here, c(t, ) and (t, ) are continuum field variables representing concentration and polarity of CSKs at time t in a two-dimensional space . F is the free energy of bundling and polarity of CSK. We employed a simple Landau energy form fc(c) = −(α2/2)(c−c0)2+(α4/4)(c−c0)4, according to which, MT accumulation is induced to form bundles. In the second term, I is the identity matrix and Q = ⊗−(1/2)Tr[⊗]I is a nematic tensor depending on the polarity of CSKs. This term represents ‘surface tension’ of the MT bundles, in which, polarity prefers parallel alignment to the CSK concentration boundary with a positive value of α, i.e., is perpendicular to the cell concentration gradient ∇c (Fig 2A(i)). The term F describes the formation of CSK polarity. By setting β(c) = β(c−c), the polarity increases with MT concentration exceeding the critical value cc, up to an amplitude of . The last term represents the tendency of the polarity to align with each other, where one constant approximation is employed. Taken together, Fcsk provides the free energy that expresses process (i), where MT concentration and polarity are coupled via nematic tensor Q. The simulation results without BBs has been shown in the methods section and S1 Fig.

F (Eq 3) and F (Eq 4) are the free energies representing processes (ii) and (iii), respectively. The variable ϕi(t,) represents the region occupied by the i-th BB, a circular disk with radius R centered at i. We set ϕi = [1+tanh(v(R−|−|))]/2 by which ϕi = 1 for the internal region of the BB, and ϕi = 0 otherwise [18]. v controls the steepness of the boundary of the BB surface and is set to be sufficiently large. F is the energy term for the excluded volume effect between BBs and CSKs, which was not considered in the previous model [3]. F is the energy term representing the excluded volume effect among BBs, to avoid overlap between different BBs, summed over all pairs of BBs.

The dynamics of the system is derived from the free energy F as follows.

In Eq (5), time evolution of CSK concentration c is derived from the energy relaxation with mass conservation, where α is the kinetic coefficient. Polymerization and depolymerization (process iv) are represented by k and k being maximum polymerization rate and depolymerization rate, respectively, with K being a constant. The mean concentration of CSKs in the system is determined by the balance between these rates. We assumed k is rescaled by k to keep the MT concentration almost constant (see Table 2). The time evolution rule of polarity, , and the center of the i-th BB, , are determined by simple relaxation of the free energy given by the derivative of F. In Eq (7), the center of the i-th BB moves by repulsive interaction with CSKs and other BBs. Stochastic noise is added to the evolution of i, where ξ is a white Gaussian random variable obtained from the normal distribution with zero-mean and standard deviation σΔt at each time step, where Δt is the discretized time step in the simulation. Finally, we assume that the polarity of CSKs guides BB orientation via the interaction between CSKs and BFs (process v). Although there is a preferred angle between the direction of CSK polarity and orientation of a BB (Fig 2A(v)), the preferred angle has been assumed to be π/2 in this model. This assumption does not change the generality of the model.

Table 2

Parameters used in Fig 2B.

Parameters
Free energy
α2	1.0	Constant of 2nd order term in Landau expansion
α4	8.0	Constant of 4th order term in Landau expansion
c0	0.5	Concentration on a steady state
D0	0.5	Diffusion constant
αQ	2.0	Anisotoropic diffusion constant along p
β	1.5	Contribution from concentration to polarity (Active term)
cc	1/3	Threshold of CSK concentration at which polarity can grow
K	3.0	Bending modulus of CSK
λcϕ	1.0	Repressive effect between BB and CSK
λϕϕ	1.0	Repressive effect between BBs
v	10.0	Steepness of the edge of BB surface
R	2.5	Radius of the disk-shaped BB
Dynamics
αJ	1.5	Constant kinetic coefficient
kp	kd(c0+KM)	Rate of polymerization of CSK*
kd	0.1	Rate of depolymerization of CSK
KM	0.2	Michaelis Menten constant in the polymerization of CSK
γP	0.5	Inverse time constant for the changing of polarity
γBB	0.3	Inverse time constant for the movement of BB
γBF	0.5	Inverse time constant for the orientation of BF
Noise strength
σBB	0.2	Standard deviation of white Gaussian random variable ξi
σBF	0.05	Standard deviation of white Gaussian random variable ξBF
System Parameters
N	80/20/30	Numbers of BBs (Fig 2, Fig 4A, 4D and 4E / Fig 3 and Fig 4B / Fig 5)
Lx	64/32/30	Length of x-axis (same as above)
Ly	64/32/45	Length of y-axis (same as above)
dx	1.0	Spacing of square grid

* k is kept to be proportional to k for the MT concentration to be almost constant.

Here, θ is the orientation of the i-th BB and = (cosθ, sinθ) is its director. The angle brackets 〈⋯〉 denote the average over the surface of the BB. Gaussian white noise term ξ is also added with standard deviation σΔt.

The above argument provides the closed form of the equations to be explored below. We have numerically solved non-dimensionalized equations where the units of time and length in the simulations correspond to 5 min and 0.04 μm, respectively. Details of the equations, parameter values, and numerical methods have been provided in the methods section.

The mathematical model reproduces BB positioning and orientation

Depending on the parameter values used in the model equations, several patterns were observed in the numerical simulations. Fig 2B shows an example of the time series obtained in the simulation when α, the coupling strength between the polarity and concentration of CSK in Eq (2), is set at α = 2.0 (see Table 2 for the other parameters). The polarity of CSK, becomes oriented in the same direction as the system and the CSK concentration, c, becomes striped along the direction of the polarity. BB arrays are formed and aligned between the stripes of CSK concentration, with all BBs being oriented in the same direction as the polarity. This pattern obtained in the numerical simulation recapitulates the experimental observations in normal conditions. We refer to this pattern as ‘alignment’. Decreasing αQ retains direction of the polarity with stripes of CSK concentration often becoming branched and BB arrays becoming winding (Fig 2C, α = 0.2), which is referred to as ‘partial alignment’ [3]. This is due to the weak coupling between polarity and CSK concentration, suggesting the importance of polarity in establishing an alignment pattern. This will be further studied in the next section.

Fig 3 illustrates the phase diagram of the patterns observed by the numerical simulation against α and the polymerization/depolymerization rate k. Some typical patterns are shown in the diagrams which are dependent on the parameter values. The boundaries between these patterns are not so sharp due to the initial condition dependence and difficulty in classification, particularly due to the appearance of mixed patterns. For weak polymerization rate (k < 10−1.5), clusters of BBs appear regardless of the value of αQ. By increasing αQ, double strand arrays of BBs are obtained. Such a pattern is occasionally observed in cultured tracheal cells [3]. On the other hand, BBs are scattered for sufficiently high polymerization rate (kd > 100.75) since the amounts of CSKs surrounding BBs recover quickly to the mean concentration and isolate them. Alignment and partial alignment patterns appear in the middle region of kd, depending on the value of αQ as mentioned above. Altogether, configuration of BBs is dependent on the polymerization rate and coupling strength between CSK polarity and concentration. Appearance of BB clusters at low polymerization rate is of interest since it resembles those observed in the nocodazole treated clusters [3], which will be discussed further.

Fig 3

Phase diagram of BB pattern against parameters αQ and kd.

αQ is the coupling constant between CSK concentration and polarity. kd is the polymerization/depolymerization rate of CSK. Patterns are classified into Alignment (A), Cluster (C), Scatter (S), partial alignment (PA), and Double strand (DS). Typical BBs pattern in each region is shown are also shown. Bars, 0.5 μm.

Polarity is required for alignment pattern of BBs

To further explore the role of polarity in CSKs, we studied the change in the behavior of the system with change in the parameters αQ and K, the latter is bending modulus of MT bundles (Eq 2). Our numerical simulation suggests that these two parameters must be sufficiently high for achieving regular alignment pattern of BBs, indicating importance of the polarity in CSKs for BB alignment pattern. In S2 Fig, the straightness of BB arrays is quantified by calculating 〈−cosψ〉, which is the average of the negative cosine of the opening angle between line segments connecting the i-th BB and its two neighboring BBs. When either αQ and K are low, arrays of BBs are short or are not aligned straight, classifying them into ‘partial alignment’ or a ‘scatter pattern’. Therefore, polarity in CSKs plays two roles in our model. It does not only convey the orientational information to BBs (Eq 8), but also enables BBs to be regularly aligned through the interaction with CSK concentration and BBs.

Bundling of CSKs enhances BB alignment formation

MTs are bundled in the apical membrane, which stiffen the apical CSK network and could contribute towards the formation of aligned BB arrays. For more insights into the role of the MT bundles for BB pattern, we changed the functional form of the free energy fc(c) in Eq (2) which drives the formation of MT bundles (condensation of MTs). We re-parametrized the free energy as fc(c) = λW(c−c0−W)2(c−c0+W)2 in which the parameter λW controls the scale of the energy (i.e., ‘depth’ of the potential), while the positions of the two minima are controlled byW. With a proper form of fc(c), BB arrays are well lined up as shown in S3 Fig, consistent with the preceding sections. Since in the system without BBs, the double minimum form is indispensable for inducing the phase separation (condensation) of CSKs resulting in the stripe pattern (see Methods and S1 Fig), we hypothesize W to be positive-finite. If that were not the case, BBs would scatter and be randomly configured. This is confirmed for large λ and small W where BBs fails to align (S3(iv) Fig). However, when both λ and W are small, there were unexpected observations in which arrays of BBs were formed and were aligned (S3(ii) Fig). Since the form of fc is flattened in the parameter region, we then conducted numerical simulation by erasing　the free energy (i.e., fc(c) = 0). We found that the aligned pattern of BBs is achieved as shown in Fig 4A. The obtained phase diagram on kd -αQ plane shown in Fig 4B is similar to those obtained in Fig 3, indicating that the condensation of CSKs is not essential for the pattern of BBs. However, the region of the ‘aligned pattern’ becomes narrower, and the BB arrays are often branched in the aligned region (Fig 4B(i)) and are less straight. As shown in Fig 4C, in the model including fc(c) term, BBs align in straighter over a wide range of noise strengths in Eq (7). Thus, MT bundles are useful for more robust and ordered positioning of BBs.

Fig 4

Model without the condensation term and simulation for nocodazole treated cells.

(A) Typical time series of BB pattern when f(c) = 0. Parameter values are the same as Fig 2B. (B) Phase diagram of the model without condensation. (i) and (ii) show examples of “alignment” and “partial alignment”, where the same parameters are used as corresponding points in Fig 3. (C) Robustness of BB for noise strength, σBB. The straightness of BB arrays, 〈−cosψ〉, was compared between models with and without f(c) term. Gray shadow area indicates σBB used in the other simulations. The error bars represent the standard error among four independent simulations. (D) Simulated patterns corresponding to nocodazole treated cells, which is produced by the model without condensation term and coupling with polarity (See Methods). Bars, 0.5 μm.

We conclude from the above numerical results that even without the condensation effect of MTs, aligned patterns can be realized to some extent by the exclusion volume effect between BBs and CSK. This was in contrast to the results of the previous model, where exclusive interaction between BBs and CSKs was not considered and condensation of CSKs was indispensable for patterning of BBs [3]. However, it is possible that the fluctuation is suppressed and the aligned pattern of BBs is achieved more robustly through the use of MT bundles. This result is biologically relevant since the maturation of multi-ciliated cells must progress under various external disturbances [10,19].

Simulation for nocodazole treated cells

Given the role of polarity and condensation of CSKs investigated in the preceding sections, here we studied the formation of BB clusters in cells treated with nocodazole, which inhibits MT polymerization [3]. As shown in Figs 3 and 4B, BB clusters are formed at low polymerization and depolymerization (k < 10−1.5), consistent with the experimental result. To clarify the mechanism of cluster formation, we simplified our model to ignore polarity and condensation of MTs, since they would play only an insignificant role when MT polymerization is inhibited. The simplified model is obtained by replacing the free energy Fcsk in Eq (2) with the following one;

The numerical simulation of the simplified model demonstrates the formation of BB clusters (Fig 4D), similar to the experimental observation in nocodazole treated cells. This result indicates that clustering of BBs occurs by merely minimizing the quadratic term of the concentration gradient, attributed to the depletion force [20]. Since the concentration gradient is localized in the boundary of BBs, accumulation of BBs is preferable for decreasing the total length of the boundary. Meanwhile, interaction length of the resulting attractive force between BBs is too short, resulting in the appearance of multiple clusters depending on the initial position of BBs.

MTs transmit directional information in PCP to BBs

All the cilia orient to the oral side of the tracheal tissue in vivo, which is coordinated by PCP. In fact, directionality of BBs is reduced by knocking out the PCP-related protein Vangl1, compared with those in wild-type (WT) [21]. Importance of PCP in cilial organization is also evident in many systems [4,6,7,22]. In Drosophila epithelial cells, the apical MTs are aligned along the proximal-distal axis (coinciding with the PCP axis), which is involved in directed molecular transportation within the cell [23-25]. Recent experiments suggest that apical CSKs mediate the transmission of directional information from PCP to BBs and cilia, where PCP proteins are localized along the oral and lung side of the adherent membrane in each cell (Fig 5A). However, the detailed information transmission mechanisms remain unclear and may be different among different cell types and species [21,26].

Fig 5

Simulation of BBs positioning by PCP.

(A) PCP proteins are localized along the oral and lung boundaries on the apical membrane. CSK concentration in the apical is low in immature cell and becomes increased during maturation. In matured cells, MTs are more distributed on the oral side (indicated by localization of DVI1/3 and Daple). (B) Polymerization rate kp is increased (blue line) and the boundary condition is changed at t = 66.7 h (gray dashed line) in the simulation. The straightness of BB arrays obtained by simulation is shown by a red line. (C) Simulation result for pattern formation of BB alignment. Black triangle represents change in boundary condition (see text). Bars, 0.5 μm.

It was observed that MT concentration is significantly high along the oral side of the cell boundary membrane in a tracheal cell, which was attributed to the interaction with PCP proteins such as Frizzled, Dishevelled, and Daple which accumulate on the oral side of the cell membrane (see Fig 5A) [3,21,27]. Although the asymmetric accumulation of the cortex MTs can also be related with MT organization over the apical membrane, and can thus proceed at the same time with the pattern formation of BBs [28], here we do not address the entire process. We restrict our study to simply demonstrate a possible mechanism of how the pattern of BBs can be modulated by the PCPs. In the framework of our model, the effect of PCP is incorporated as the boundary condition as follows [3]. For the oral side boundary, concentration of CSKs is set to be high (i.e., Dirichlet condition c = cH is applied, see Methods for details) and direction of the polarity is restricted to be parallel along the border. Neumann conditions are applied for the variables for the other boundary. In order to demonstrate the maturation process of a multi-ciliated cell during high CSK concentration, we simulated a condition where the CSK concentration on the oral boundary is switched from Neuman to Dirichlet boundary conditions mentioned above, at time t = 66.7 h, and the parameter kp also increases with time (Fig 5B blue line; see Methods for details).

The time series of the simulation result is shown in Fig 5C. Initially scattered pattern of BBs (t = 10 h in the panel) evolve to the partial alignment pattern (t = 50 h and 80 h), where short or branched arrays of BBs are formed by the increase in CSK concentration and BBs orient to the same direction locally (Fig 5C). BBs eventually establish a fully aligned pattern that is parallel to the oral border through connection and reconnection of BB arrays, (t = 110h). This time evolution of BB pattern agrees with those observed in the normal maturation process of a multi-ciliated cell.

Discussion

On the apical membrane of epithelial cells, CSKs such as MTs, actins, and 10 nm-filaments can perform various biological functions through the process of pattern formation attributed to the self-organizational ability of CSKs [15,17,29-31]. In the present study, we investigated positional and orientation orders of BBs in the apical membrane of multi-ciliated cell which are responsible for coordinated ciliary beating to generate mucociliary transport. Apical MTs, which are bound to the BF of BBs, play a crucial role in achieving the positional and orientational orders of BBs. Most importantly, consideration of the symmetry led us to propose that polarity should be maintained in the MT bundle filaments for properly guiding orientation in BBs. This idea is supported by the biased distribution of angles between BB orientation and MT filaments in TEM images (Fig 1D). Although more direct evidence of the presence of polarity in MT bundles is still lacking, this can be investigated and confirmed in future experiments.

We propose a mathematical model for patterning of BBs, taking into consideration the polarity of MT bundles. The model was constructed based on assumptions such as CSK characteristics and volume exclusion among BBs and CSKs. Since the fundamental mechanism of how BBs align and orient is not yet known, the model was not aimed at being precise or quantitative, but rather provided a basic description of the processes for understanding the underlying mechanism following the generic formalism [13-15]. The derived model is distinct from the previous one [3] in incorporating polarity as well as the repulsive interaction between BBs and CSKs, by which it can reproduce both the positional and orientational orders of BBs in multi-ciliated cells. Analysis of the model demonstrated that polarity is crucial not only for determination of BB orientation, but also for achieving alignment pattern of BBs, the latter depending on the nematicity of the CSK determined by polarity. It was also shown that BBs can align even without the condensation effect of CSK concentration, i.e., exclusion of BBs and CSKs was sufficient for the formation of alignment pattern. However, the formation of MT bundles enabled robust formation of the alignment pattern against external disturbance. These simulation results are feasible given the experimental observation that MTs are indispensable for establishing alignment pattern in the apical membrane of multi-ciliated cells where MT bundles segment the aligned array of BBs. Positional order of BBs was disrupted by suppressing MT polymerization with BB clusters being formed, which could also be reproduced by the present model. Taken together, it was demonstrated that incorporation of CSK polarity in the model was necessary and sufficient to reproduce prior experimental observations.

It should be noted that our model still depends on various assumptions and misses some crucial aspects observed in experiments. First, in our model, the effect of PCP was simplified into a fixed boundary condition (i.e., Dirichlet condition) that represents accumulated tyrosinated MTs on the oral side of cellular cortex [21]. As explained before, this accumulation of the cortical MTs itself can be regulated by both PCP proteins and apical MTs. These processes should be incorporated to study a more precise time course of patterning of BBs. Second, in our model, it was assumed that the orientation of BBs merely follows the direction of MTs (process (v) in the model assumption) for which the detailed mechanism was not taken into account. Since BF contains γ-tubulin [9], it will be interesting to explore whether nucleation of the apical MTs at BF has a possible impact on the determination of BB orientation to MTs. In addition, while PCP and apical CSKs are indispensable for the coordination of BBs and ciliary beating, it was reported that hydrodynamic interaction or direct physical contacts among cilia were part of a positive feedback loop for accurate BB patterning and contributed to the correct directionality in BBs and ciliary beating to enhance mucociliary transport [19,26,32]. It would be interesting to extend our model to understand how the mechanical perturbation from the synchronously beating cilia to the BBs in apical membranes enables smooth BB rearrangement and promotes their alignment pattern.

Finally, we would like to discuss the implication of our study in determining cellular chirality. Presence of polarity in CSKs aligning perpendicular to the PCP axis (oral to lung axis) indicates that the cell can distinguish between left and right directions from the PCP axis (Figs 1A and 5A, blue arrows). From a theoretical point of view, our study provides a possible mechanism for implementing cellular chirality by using PCP and CSK polarity. In the biological system, chirality appears at various levels ranging from molecular [33] and cellular [34] to the tissue level [35]. Cellular chirality is responsible for body formation. For example, cell intrinsic chirality underlies left-right asymmetry of chicken cardiac looping [36,37]. Nodal cilia are well-known examples in which the cilia are tilted in a definite direction and generate the left-ward flow on the nodal tissue surface, which induces the left-right asymmetry of animal body [7,38,39]. Although the mouse trachea studied in the present work has no apparent chiral asymmetry, it is possible that the system can regulate cellular chirality through some mechanisms, which is of potential importance in generating chirality at the tissue level. We believe that our model provides useful insights into exploring the common mechanisms of cellular and tissue chirality.

Methods

Mathematical model for BB pattern formation in multi-ciliated cells

Details of the model Eqs (5–7) in the main text are as follows.

Ginzburg-Landau potential has been employed for this model which has two minima at . The last term in Eq (11) represents the reactive term to align the polarity along the CSK concentration gradient, where M = (∇c)⊗(∇c)−(1/2)Tr[(∇c)⊗(∇c)]I. The reactive term associated with β(c) was ignored for simplicity in Eq (10). This was verified by interpreting that the emergence of polarity is governed by active dynamics. In practice, the results are almost independent of the presence of the reactive term. In Eq (3), we adopted F = λϕc2 for exclusive volume effect between i-th BB and MTs. Alternative symmetrically permissible form is to adopt , by which we found similar patterns with appropriate choice of and λ.

The dynamics of CSK and BBs was numerically simulated by integrating the differential Eqs (8 and 10–12). The non-dimensionalized equations of CSK defined on a rectangular domain Lx × Ly were calculated. The time integration is done by the fourth-order finite difference scheme with time step Δt = 1.0 × 10−3. Periodic boundary condition was employed unless stated otherwise. The system size was chosen as Lx = Ly = 64, for which the domain is discretized by a square grid with Δx = 1.0 spacing. The number of BBs is N = 80 for most of the simulations. Numerical simulations were performed with system size Lx = Ly = 32 and N = 20 for producing the phase diagrams in Figs 3 and 4. Four independent simulations were performed for each set of parameters to check the dependency on initial conditions. Initial conditions were set as follows. The position and orientation of BBs, i and θ, were random, the polarity of CSK was set at || = 0.1 with random direction, and the concentration of CSK is set at c = 0.2 with small value of additional noise.

The values of the parameters used in the simulations are summarized in Table 2. Since exact kinetic values were not determined experimentally, we determined them so as to produce experimental patterns, using the following rationale. The units of length and time were set to be equal to 0.04 μm and 5 min, respectively. Radius of a BB was chosen as R = 2.5 corresponding to 0.1 μm [3,11]. System size was Lx = Ly = 64 × 64 corresponding to 2.56 × 2.56 μm2, approximately quarter of the apical area of single multi-ciliated cells. N = 80 since ~300 BBs exist per cell. Typical speed of BBs was about 0.02–0.03μm/ 5min [3]. Hence, γ was chosen to satisfy dX/dt~γR = 0.03μm/ 5min. The concentration of CSKs was not known. Therefore, we normalized it and used the dimensionless concentration for c by setting parameters to ensure c- = 0 and c+ to be almost unity. Amplitude of polarity || was also chosen to be almost unity. Another prerequisite for c, modeled using Ginzburg-Landau potential combined with the reaction term [40], was that the system must show bundling. Reaction rate kd should be smaller than the critical value evaluated by linear stability analysis, (see S1 Fig). Typically, we used which corresponds to ~15 min in reality. This is comparable with the observed time scale of MT bundle dynamics in vitro [41]. In addition, we set α = γK following the model for actin filament in [15]. With these choice of parameters, the pattern formation of BBs took 600~1200 time in the simulations (Figs 2 and 5), corresponding to 2~4 days. This is consistent with actual developing time and experimental observation in cultured tracheal cells [3]. We should note that the parameter values themselves, such as equilibrium concentration of MTs (kp and kd), can be gradually changed during the developmental stage, which could govern the time scale of patterning ultimately.

We set kd = 1.0 × 10−2 to obtain the results in Fig 4D, where polarity is ignored and the dynamics of CSK and BBs are simulated by using Eq (12) and

We set the system size as N = 30, Lx = 30, and Ly = 45 to simulate the model coupled with PCP information shown in Fig 5. As briefly mentioned in the main text, Neumann boundary condition was initially employed, followed by time step t = 66.7 h. The boundary condition of the oral side (top boundary in Fig 5C) was changed to Dirichlet where CSK concentration was fixed at c = cH to mimic the experimental observation that the distribution of MTs is dense in the opposite side of the cell boundary where PCP protein Vangl1 is accumulated (i.e., lung side) [3,21]. cH was chosen as . Polymerization rate kp was set to be an increasing function of time as kp/k = 0.6 + 0.1 tanh(1.0 × 10−3t−0.5) (Fig 5B, blue line), since MTs and intermediate filaments increase during maturation of cells [12].

Model without BBs

It is useful to summarize the case where BBs are absent, to understand the behavior of the model. The model is composed of phase-separation dynamics with first order chemical reaction [42] and polarity dynamics [13]. The phase separation is driven by Landau-Ginzaburg potential f(c). We have numerically solved Eqs (10 and 11) without the terms representing interaction with the BBs.

As shown in S1 Fig, the pattern of CSK concentration is dependent on αQ and the polymerization/depolymerization parameter kd [42]. For lower chemical reaction rate (i.e., smaller kd), the CSK domain concentration is larger. CSK concentration shows stripe pattern in the middle range of kd. The range of aligned stripe is wider when the coupling constant with polarity, αQ, is larger, while it is very narrow at αQ = 0. In addition, the stripe pattern is less sharp for smaller αQ. The CSK concentration becomes uniform for larger value of kd. In summary, the model shows aligned robust bundles of CSKs for proper value of kd and higher value of αQ. The model without BBs satisfies the minimal prerequisite to express the bundling of CSKs in our model.

Patterns of CSKs without BBs.

The patterns are dependent on αQ and kd; Striped pattern appears for finite value of αQ (striped patterns which are not sharp may appear at αQ = 0 in a narrow parameter region). CSKs concentration becomes uniform at very high values of polymerization/depolymerization rate kd. Bars, 0.5 μm.

(TIF)

Click here for additional data file.

BB positional order depends on the elasticity parameter, K.

Color indicates the straightness of BB quantified by 〈−cosψ〉. Inset figures show typical pattern of BBs. Gray dashed line (K = 3.0) corresponds to the line connecting between symbol A and PA in Fig 3. Bars, 0.5 μm.

(TIF)

Click here for additional data file.

Dependence of BB pattern on the form of condensation term fc(c).

The form of fc is expressed as fc(c) = λW(c−c0−W)2(c−c0+W)2. Left: phase diagram against λW and W. Color indicates the straightness of BB (〈−cosψi〉). Right: Four representative patterns are shown (i-iv) with corresponding black squares in the left panel. Bars, 0.5 μm.

(TIF)

Click here for additional data file.

Model simulation of BBs corresponding to Fig 2B.

BBs reach “alignment” pattern. Bars, 0.5 μm.

(MP4)

Click here for additional data file.

Model simulation of BBs corresponding to Fig 2C.

BBs reach “partial alignment” pattern. Bars, 0.5 μm.

(MP4)

Click here for additional data file.

Model simulation of BBs without the condensation term, corresponding to Fig 4A.

Bars, 0.5 μm.

(MP4)

Click here for additional data file.

Model simulation for nocodazole treated cells, corresponding to Fig 4D.

Bars, 0.5 μm.

(MP4)

Click here for additional data file.

Model simulation of BBs aligned by PCP, corresponding to Fig 5C.

Bars, 0.5 μm.

(MP4)

Click here for additional data file.

7 Sep 2019

Dear Dr Ishihara,

Thank you very much for submitting your manuscript 'Cytoskeleton polarity is essential in determining direction in basal bodies of multi-ciliated cells' for review by PLOS Computational Biology. Your manuscript has been fully evaluated by the PLOS Computational Biology editorial team and in this case also by independent peer reviewers. The reviewers appreciated the attention to an important problem, but raised some substantial concerns about the manuscript as it currently stands.

In particular, the model presented is a next step over your previous model presented in Ref. [3], in that the new model considers the polarity of the CSKs. Your new work shows that if the polarity of the CSKs, i.c. microtubules (MTs) is considered, this can explain how the BB arrays become orientationally ordered. However, reviewer 1 notes that such polarity/parallel ordering is likely not found - the MTs are ordered antiparallel with a slight bias at best. The reviewer argues that, in cases where such slight MT orientation is found, it is likely to be the result, not the cause of BB orientation. Also if cilia motion is disrupted, thus disrupting fluid flow, the orientation of the BBs is lost; although this key role of fluid flow is mentioned in the discussion, the model does not provide an explanation for this discrepancy with the experimental observations.

While your manuscript cannot be accepted in its present form, we are willing to consider a revised version in which these issues and the other issues raised by the reviewers have been adequately addressed. We cannot, of course, promise publication at that time.

Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out.

Your revisions should address the specific points made by each reviewer. Please return the revised version within the next 60 days. If you anticipate any delay in its return, we ask that you let us know the expected resubmission date by email at ploscompbiol@plos.org. Revised manuscripts received beyond 60 days may require evaluation and peer review similar to that applied to newly submitted manuscripts.

In addition, when you are ready to resubmit, please be prepared to provide the following:

(1) A detailed list of your responses to the review comments and the changes you have made in the manuscript. We require a file of this nature before your manuscript is passed back to the editors.

(2) A copy of your manuscript with the changes highlighted (encouraged). We encourage authors, if possible to show clearly where changes have been made to their manuscript e.g. by highlighting text.

(3) A striking still image to accompany your article (optional). If the image is judged to be suitable by the editors, it may be featured on our website and might be chosen as the issue image for that month. These square, high-quality images should be accompanied by a short caption. Please note as well that there should be no copyright restrictions on the use of the image, so that it can be published under the Open-Access license and be subject only to appropriate attribution.

Before you resubmit your manuscript, please consult our Submission Checklist to ensure your manuscript is formatted correctly for PLOS Computational Biology: http://www.ploscompbiol.org/static/checklist.action. Some key points to remember are:

- Figures uploaded separately as TIFF or EPS files (if you wish, your figures may remain in your main manuscript file in addition).

- Supporting Information uploaded as separate files, titled Dataset, Figure, Table, Text, Protocol, Audio, or Video.

- Funding information in the 'Financial Disclosure' box in the online system.

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, https://pacev2.apexcovantage.com PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at figures@plos.org.

To enhance the reproducibility of your results, we recommend that you deposit your laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. For instructions see here.

We are sorry that we cannot be more positive about your manuscript at this stage, but if you have any concerns or questions, please do not hesitate to contact us.

Sincerely,

Roeland M.H. Merks, Ph.D

Associate Editor

PLOS Computational Biology

Mark Alber

Deputy Editor

PLOS Computational Biology

A link appears below if there are any accompanying review attachments. If you believe any reviews to be missing, please contact ploscompbiol@plos.org immediately:

[LINK]

Reviewer's Responses to Questions

Comments to the Authors:

Please note here if the review is uploaded as an attachment.

Reviewer #1: The manuscript by Namba and Ishihara deals with the fascinating question of how cilia in a multiciliary array become coordinately polarized. The authors provide a mathematical model for how this is achieved that relies heavily on the interactions between basal bodies and the microtubule cytoskeleton. The authors do an admirable job of summing up the relevant literature and there is certainly some merit to their model. However, I have several serious concerns that warrant consideration.

1. MT bundles can be parallel or anti parallel. The information in the literature would suggests that these MT bundles are mixed or anti parallel, but the current equations if I understood them correctly rely on parallel bundles all going in the same direction. This is not likely to be the case. At the very best there will be a slight bias in the directionality of the MTs. I think this mixed array needs to feature more heavily into any equations.

2. Similar to above the authors state “Although more direct evidence of the presence of polarity in MT bundles is still lacking, this can be investigated and confirmed in future experiments.” Multiple groups have tried to do this in multiple systems. There is simply not a strong directional bias to the MT orientation before BB polarization. Vladar et al. 2013, Boutin et al 2014, Kim et al 2018 all show some polarized accumulation of MTs in MCCs but this comes after BB polarization and likely reflects the bias in orientation of MTs nucleating from the BF that are all pointing in the same direction. If oriented MTs were defining the polarity for BB then it would come first, which does not seem to be the case. Additionally, all of the mentioned data have polarized subset of MTs parallel with the axis of orientation (and primarily at the cortext) whereas this model is based on the pool of MTs interacting with BBs which must be more perpendicular to the axis of orientation.

3. The authors state “However, this binding between BF and MTs allows BBs to associate with either of the two neighboring MT bundles and hence they cannot distinguish between the two opposite directions (Fig 1C (i)). Therefore, how BBs can determine their definite direction remains a mystery.” If I understand things correctly, the authors are treating the MTs and the BB independently. Figure 1C suggests that the BF are interacting with the MTs, which is somewhat consistent with the EM literature. However, what I think is missing from this discussion is where the MTs are originating. I think there is strong indications in the literature that the MTs are in fact at least in part nucleating from the BF. Therefore it is not really appropriate to consider these things completely independent of one another.

4. Related to above all MTs in this manuscript are created equal which is almost certainly a critical oversimplification. The evidence from Vladar et al. would suggest that there is a distinct pool of MTs labeled with tyrosinated tubulin that interact with cell cortex in a polarized manner. However the majority of MTs are not tyrosynated in these cells.

5. Finally and perhaps most concerning is the lack of integration of the key information that hydrodynamic forces are critical to the polarization of these BBs. While the authors mention this in their discussion it is completely lacking from the model. We know from multiple systems that if cilia motility is disrupted (with the rest of the MT cytoskeleton intact) then BBs will not orient. This tells us that the MT BB interactions alone are not sufficient to drive orientation. While I appreciate that this is a complicated question and the connection of MTs and BBs is an important aspect, a model that does not explain the basic experimental evidence needs improvement.

Minor:

“In particular, the attachment of BF to MTs might be important for the determination of BB direction, but its mechanistic details have not been identified as yet.” I more or less agree with this statement and also feel that in general the authors have done a thorough job of referencing the literature, however I think referencing the Turk et al paper on gamma tubulin at the basal foot might be appropriate as this is probably the closest we have seen to mechanism for nucleating MTs from the BF.

Reviewer #2: Uploaded attachment

**********

Have all data underlying the figures and results presented in the manuscript been provided?

Large-scale datasets should be made available via a public repository as described in the PLOS Computational Biology data availability policy, and numerical data that underlies graphs or summary statistics should be provided in spreadsheet form as supporting information.

Reviewer #1: Yes

Reviewer #2: Yes

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

Submitted filename: PLOScomments.pdf

Click here for additional data file.

28 Nov 2019

Submitted filename: Response_to_reviewers.pdf

Click here for additional data file.

9 Jan 2020

Dear Dr Ishihara,

We are pleased to inform you that your manuscript 'Cytoskeleton polarity is essential in determining orientational order in basal bodies of multi-ciliated cells' has been provisionally accepted for publication in PLOS Computational Biology.

Before your manuscript can be formally accepted you will need to complete some formatting changes, which you will receive in a follow up email. Please be aware that it may take several days for you to receive this email; during this time no action is required by you. Once you have received these formatting requests, please note that your manuscript will not be scheduled for publication until you have made the required changes.

In the meantime, please log into Editorial Manager at https://www.editorialmanager.com/pcompbiol/, click the "Update My Information" link at the top of the page, and update your user information to ensure an efficient production and billing process.

One of the goals of PLOS is to make science accessible to educators and the public. PLOS staff issue occasional press releases and make early versions of PLOS Computational Biology articles available to science writers and journalists. PLOS staff also collaborate with Communication and Public Information Offices and would be happy to work with the relevant people at your institution or funding agency. If your institution or funding agency is interested in promoting your findings, please ask them to coordinate their releases with PLOS (contact ploscompbiol@plos.org).

Thank you again for supporting Open Access publishing. We look forward to publishing your paper in PLOS Computational Biology.

Sincerely,

Roeland M.H. Merks, Ph.D

Associate Editor

PLOS Computational Biology

Mark Alber

Deputy Editor

PLOS Computational Biology

Reviewer's Responses to Questions

Comments to the Authors:

Please note here if the review is uploaded as an attachment.

Reviewer #1: In general I think modeling studies are really important and so I like to be supportive. However, I also get very frustrated by them as they are always over simplifications. I think this study is reasonable for what it is but they have not been able to deal with my primary concern that their model does not incorporate what the field knows is a critical aspect, namely hydrodynamic forces. I appreciate that this would be challenging and they certainly mention it and propose to add it in future. But as a researcher in the field I am not sure this work changes anything about the way I think about this question. The discussion on chirality is interesting but seems out of place in the current work. Overall, I think this work does a good job of addressing some of the data in the field and thus it does have merit, and it is clear that it is hard to isolate the importance of individual factors when we know that these are part of a larger combined effort. I think the authors should be commended on their effort and I think it is a worthwhile contribution that may be the foundation for future more inclusive models. As such I think it has value to the field and I am supportive of publication.

Reviewer #2: The manuscript is considerably improved in terms of readability. I enjoyed reading it a second time. The manuscript is well and very carefully written, including multiple justifications and considerations that can serve as starting work for future studies. The introduction and discussions have also been enriched. I'm happy to recommend it for publication.

**********

Have all data underlying the figures and results presented in the manuscript been provided?

Large-scale datasets should be made available via a public repository as described in the PLOS Computational Biology data availability policy, and numerical data that underlies graphs or summary statistics should be provided in spreadsheet form as supporting information.

Reviewer #1: Yes

Reviewer #2: Yes

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

13 Feb 2020

PCOMPBIOL-D-19-01166R1

Cytoskeleton polarity is essential in determining orientational order in basal bodies of multi-ciliated cells

Dear Dr Ishihara,

I am pleased to inform you that your manuscript has been formally accepted for publication in PLOS Computational Biology. Your manuscript is now with our production department and you will be notified of the publication date in due course.

The corresponding author will soon be receiving a typeset proof for review, to ensure errors have not been introduced during production. Please review the PDF proof of your manuscript carefully, as this is the last chance to correct any errors. Please note that major changes, or those which affect the scientific understanding of the work, will likely cause delays to the publication date of your manuscript.

Soon after your final files are uploaded, unless you have opted out, the early version of your manuscript will be published online. The date of the early version will be your article's publication date. The final article will be published to the same URL, and all versions of the paper will be accessible to readers.

Thank you again for supporting PLOS Computational Biology and open-access publishing. We are looking forward to publishing your work!

With kind regards,

Sarah Hammond

PLOS Computational Biology | Carlyle House, Carlyle Road, Cambridge CB4 3DN | United Kingdom ploscompbiol@plos.org | Phone +44 (0) 1223-442824 | ploscompbiol.org | @PLOSCompBiol