Predictive model of spatial nematic order in confined cell populations

Introduction

In the fields of tissue engineering and biotechnology, fabrication technologies for two-dimensional cellular tissues, including cell sheet technology, have been developed for the quantitative analysis of tissues, disease modeling, drug discovery, and regenerative medicine technology^1,2,3. Much of the research has focused on cell alignment in two-dimensional cellular tissues because many kinds of spindle-shaped cells, such as fibroblasts, myocytes, and epithelial cells, show alignment order in two-dimensional tissues when they are in a confluent state. In addition, cell alignment is related to many biological phenomena; for example, it is known that well-aligned myoblast cells enhance their differentiation efficiency from myoblasts to myotubes^4,5,6 and that aligned cardiomyocyte cells can generate larger contraction, which is an important factor for fabricating engineered muscular tissues^7,8.

To control cell alignment, many researchers have focused on topographical guides^9,10,11,12 since adhesive cells tend to align along the geometrical shape of topographical guides. However, this approach has a problem in that topological defects, at which cell alignment cannot be defined, are spontaneously generated under some conditions, which will lead to a disturbance of the cell alignment. Duclos et al. first reported that two +1/2 topological defects, around which the cell alignment rotates +180°, are generated in circular domains¹³, and Endresen et al. reported defects of integer charges, around which the cell alignment changed ±360°, were induced by topographical patterns¹⁴. Some research has found that topological defects are more localized in non-circular domains due to their asymmetry^15,16. The topological defects can also cause mechanically unstable states in cellular tissues; for example, Kawaguchi et al. revealed that neural stem cells tend to aggregate around defects¹⁷, and Saw et al. showed that topological defects could induce apoptosis in Madin–Darby canine kidney cells (MDCK)¹⁵. Guillamat et al. reported myoblast cells confined in a small circle generated full-integer defects, which triggered differentiation into myotubes¹⁸. These studies have demonstrated that the spatial distribution of topological defects is an important physical factor for cell tissue uniformity. In addition, topographical control of the defects has also been applied to the estimation of mechanical parameters, such as polar order parameter, flow alignment parameter, and rotational viscosity of cellular tissues, by measuring cell alignment and velocity around defects¹⁹.

To predict the behavior of topological defects of cellular tissues in confined geometries, it is necessary to model and calculate the localization of the defects. One promising approach for such calculations is modeling cells as nematic liquid crystals^{20,21,22,23,24}. As Duclos et al. discussed¹³, spindle-shaped cells show a tendency of self-alignment as nematic liquid crystals, and the nematic alignment order becomes high as the cell density becomes high. During this process, the +1/2 defects in the cell alignment behave like self-propelled particles that decrease the Frank elastic energy of the cell monolayer, which has been experimentally and theoretically confirmed for the circular domains¹³. To establish a theoretical model of the defect positions, Duclos et al.¹³ derived a closed formula to calculate the elastic energy of nematic liquid crystals confined in circular domains. By minimizing the elastic energy with respect to a defect position, they showed that the elastic energy is minimized when the defects are symmetrically located with respect to the center of the circle and the radial position is 5^−1/4R, where R is the radius of the circle. This theoretically predicted value was consistent with the mean value of the defect positions in the experiments.

As described above, the movement of the defects is driven by the activity of the cell alignment, which is not considered in non-living liquid crystals. Since the activity of the cell alignment is not deterministic, we have to consider the randomness of the defect movement to quantify the spatial distributions of defect position in experiments. To this end, Duclos et al. also introduced a probability density function defined from the elastic energy and the effective temperature¹³. The effective temperature expresses the magnitude of the fluctuations of the defect positions. They estimated the effective temperature from the experimental data and found that it was relatively low, which implied the activeness of the defect motion was damped by friction between the cells and substrate. Thus, the spatial order of the cell alignment and defects in circular domains can be quantitatively analyzed and predicted by an explicit calculation of the elastic energy and the probability density of the defects. The estimated effective temperature is especially important for quantifying how much the defect motion is damped and predicting the resultant spatial distribution of the defects. Considering that existing studies showed that defects are localized in asymmetric domains^15,16, it is desirable to generalize this calculation approach to arbitrary-shaped geometries, which will allow the quantitative analysis of defect distributions for asymmetric domains and predictive design of geometries for controlling defect positions in preparation for experiments. However, the existing study¹³ achieved explicit calculations of the elastic energy and probability density function only for circular domains and made a strong assumption that the defects are located at the same radial positions, which cannot be assumed for arbitrary-shaped domains.

In this paper, we propose a new calculation method to quantify and predict cell alignment and spatial distributions of defect positions in elongated cell populations confined to arbitrary, simply connected domains by combining Duclos et al.’s modeling of defect distributions for circular domains¹³ with our previous study on explicit formulas for cell alignment in arbitrary-shaped simply connected domains²⁵. To this end, we extend our previous theoretical work²⁵, which proposed an explicit formula to calculate the alignment angle and elastic energy of nematic ordered cells in simply connected domains based on complex potentials and conformal mappings. First, we present experimental results for defect localization in several non-circular domains parametrized by conformal mappings. Then, we introduce a stochastic differential equation to model the fluctuation of the defect positions based on previous research¹³ and propose a method to explicitly express the equilibrium distribution. Finally, we propose a method to estimate the parameters of the stochastic fluctuations of the defect positions based on a Markov chain Monte Carlo (MCMC) method and verify the proposed method by numerical simulations and culture experiments.

Results

Topological defects in confined cell populations on microwells

We fabricated microwell arrays from polydimethylsiloxane (PDMS) (Fig. 1(a)) for culturing adhesive cells in confined geometries and forcing the cells to align along the wells. To simplify the following theoretical investigation, we designed the geometries of the microwells by using conformal mapping with a single parameter. Specifically, the domain of the well was parametrically expressed as the image of a unit disk ({mathbb{D}}:={zin {mathbb{C}},| ,| z{| }^{2} < 1}) under a conformal mapping w = g(z) = z + αz³, where 0 < α < 1/3. The shape of the domain became constricted in the vertical axis (Fig. 1a). The depth of the microwells was 45 μm and the area of the domain was scaled to be equal to that of a circle with a radius of 300 μm. We seeded and cultured mouse myoblast (C2C12) cells on the microwells (see the Methods section for details). As the cells divided and migrated on the wells, the cells tended to be elongated along the boundaries of the wells. Finally, the elongated cells completely covered the microwell and generated two +1/2 defects, around which the cell orientation changed +180° in a counter-clockwise manner (Fig. 1b).

Predictive model of spatial nematic order in confined cell populations — **Fig. 1: Cell alignment and defect distributions for C2C12 cells on PDMS substrates.**

As described in the Methods section, we observed the cell alignment as phase-contrast images 1 day after changing the growth medium to the differentiation medium, which meant that the cell division was suppressed and the cell alignment was almost constant. By analyzing the cell orientation for each well, we detected the defect positions in each well and then acquired spatial distributions of the defects in an almost equilibrium condition by merging the detected defect positions (blue points in Fig. 1c). Figure 1c shows that the defect positions tended to be non-uniformly scattered around certain points, which are the minimizers of the Frank elastic energy, explained in a subsequent section (red points in Fig. 1c). This localization depended on the parameter α and became more focal as α increased, which implied that the microwell shapes affected the localization of the defect positions.

To quantify the average behavior and the spatial fluctuations of the defect positions, we calculated the mean and standard deviation of the defect positions from the experimental results as the statistics of the spatial distributions of the defect positions in the equilibrium state, which were then compared with the numerically optimized defect positions. For the defects in the right half-plane, the horizontal (u-axial) and vertical (v-axial) components of their positions were defined as u(>0) and (utan phi), respectively, where ϕ is the angle of the line between the origin and the defect measured counterclockwise from the positive horizontal axis (Fig. 2a, b). Considering the line symmetry of the geometries with respect to the v-axis, the defects in the left half-plane were reflected with respect to the v-axis, and then the components of their positions were calculated. As a result, the mean value of the horizontal component u tended to be larger as α increased (blue line in Fig. 2c) and its standard deviation tended to decrease (blue line in Fig. 2e). In contrast, the mean value of the angle ϕ was almost zero for all the parameters (blue line in Fig. 2d), but its standard deviation was still decreasing with increasing α (blue line in Fig. 2f), which implied that the defects were localized around the horizontal axis (u-axis) rather than the vertical axis (v-axis). The numerically predicted defect positions, as explained in the next section, also agreed with the trend of the mean values with respect to the increase of α (red lines in Fig. 2c, d).

**Fig. 2: Statistical analysis of the defect positions from the experimental results.**

As shown in Fig. 1, the two + 1/2 defects were localized around the minimizer of the Frank elastic energy. Compared to Duclos et al.¹³, which considered circular domains, defects were anisotropically localized in our experiments as α increased, which meant that the geometry became non-circular and constricted. The average trend of the defect positions in the experiments was consistent with the numerically predicted defect positions based on our previous study²⁵. However, our previous method could not evaluate the degree of the localization around the minimizer of the Frank elastic energy, because the previous method could only calculate the elastic energy. To explain why the defects were more localized as α increased, the following sections are devoted to a detailed analysis of the defect distributions by building a stochastic model of defect positions.

Physical model of topological defects based on nematic liquid crystal theories

In this section, we review the physical model of elongated cells as nematic liquid crystals following the previous study²⁵, used as a basis for the stochastic model in the next section. We take Ω as the microwell domain, expressed as the image of a unit disk ({mathbb{D}}) under the conformal mapping g(z) = z + αz³. We consider two +1/2 defects in Ω with positions w_k = g(z_k) (k = 1, 2), where ({z}_{k}in {mathbb{D}}) is the corresponding point of the kth defect in ({mathbb{D}}). Under a boundary condition that the cells align along the boundary of the domain ∂Ω, the cell angle at (u, v) ∈ Ω can be calculated by taking the real part of the following complex potentials f_Ω(w ∣ w₁, w₂) according to our previous study²⁵:

$${f}_{Omega }(w,| ,{w}_{1},{w}_{2})={f}_{{mathbb{D}}}({g}^{-1}(w),| ,{z}_{1},{z}_{2})+h({g}^{-1}(w)),$$

(1)

where

$${f}_{{mathbb{D}}}(z,| ,{z}_{1},{z}_{2})=-frac{i}{2}mathop{sum }limits_{k=1}^{2}(log (z-{z}_{k})+log (1-overline{{z}_{k}}z))-frac{pi }{2},$$

(2)

and

$$h(z)=-ilog g^{prime} (z)=-ilog (1+3alpha {z}^{2}),$$

(3)

where (overline{{z}_{k}}) is the complex conjugate of z_k. The first term ({f}_{{mathbb{D}}}(z,| ,{z}_{1},{z}_{2})) in Equation (1) represents the cell alignment in the unit disk ({mathbb{D}}) and the second term h(z) is the complementary term due to the conformal mapping g(z). It is noted that h(z) was expressed based on the Poisson–Schwartz integral formula in the previous study²⁵, while we found that h(z) can be explicitly represented by the conformal mapping g(z). The derivation of the holomorphic function h(z) is described in the Methods section.

Figure 3 shows an example calculation of the cell alignment for α = 0.2. The complex potential h(z) can be interpreted as the sum of the two potentials corresponding to the two topological defects with topological charge +1, which are located at (z=pm frac{1}{sqrt{3alpha }}i). Thus, the cell alignment in non-circular geometries is determined by mirror images of the original defects and +1 defects located at zero points of (g^{prime} (z)) (Fig. 3b). In general, a pair of defects is subject to repulsive or attractive forces when the topological charges of the two defects have the same or opposite signs, respectively. Since the defects in the unit disk are subject to repulsive forces caused by both the +1 defects fixed on the imaginary axis and the mirror images of the original defects, we can find the optimal defect positions which balance these repulsive forces. According to our previous results²⁵, the elastic energy stored by these repulsive forces, known as Frank elastic energy under one-constant approximation²⁶, can be calculated from the following contour integral:

$$F=frac{K}{2}{iint}_{tilde{!Omega }}{left| frac{partial {f}_{Omega }(w)}{partial w}right| }^{2}dudv$$

(4)

$$=frac{K}{2}cdot frac{1}{2i}{oint }_{!partial {g}^{-1}(tilde{Omega })}frac{partial {f}_{Omega }(g(z))}{partial z}overline{{f}_{Omega }(g(z))}dz,$$

(5)

where K is Frank’s elastic constant, (tilde{Omega }:=Omega backslash ({D}_{1}(epsilon )cup {D}_{2}(epsilon ))), D_k = {w_k + ϵe^iθ ∣ θ ∈ [0, 2π)}, and ϵ is a small positive constant which represents the size of the point defect. In Equation (5), (overline{{f}_{Omega }(g(z))}) is the complex conjugate of f_Ω(g(z)). Then, the repulsive force acting on the kth defect can be represented as (-frac{partial F}{partial overline{{z}_{k}}}) in the unit disk ({mathbb{D}}). Since the repulsive forces are balanced when the elastic energy F is minimized with respect to the defect positions, it is expected that the defects are generated at the minimizers of F. By assuming that the defects are placed at point-symmetrical positions with respect to the origin and calculating the elastic energy for various defect positions, the minimizers for each well geometry can be calculated as shown by the red dots and lines in Figs. 1c and 2a, b, respectively. The u-component of the minimizer becomes larger as α increases, and the minimizer is located on the u-axis, which is consistent with the experimental results.

**Fig. 3: Numerical example of calculating the cell alignment in the domain expressed by the conformal mapping g(z) = z + αz³.**

Stochastic model of defect positions

In this section, we propose an explicit expression of the spatial distribution of the defects to explain the localization of the defects shown in Figs. 1, 2. For this, we follow Duclos et al.’s study¹³ and model the stochastic fluctuations of the defects as random walks in an external force field, since the defects at w = w_k = u_k + iv_k (k = 1, 2) are considered to be subject to the force (-frac{partial F}{partial overline{{w}_{k}}}). Let ({boldsymbol{w}}:={[{u}_{1},{v}_{1},{u}_{2},{v}_{2}]}^{T}) be a vector of the defect positions. We also denote a multivariate Gaussian process as R(t): = [R₁, R₂, R₃, R₄]^T, where each component R_k (k = 1, ⋯, 4) is independent of the other components and satisfies the following conditions:

$$begin{array}{ll}{mathbb{E}}[{R}_{k}],=,0,\ {mathbb{E}}[{R}_{k}(t){R}_{k}(t^{prime} )],=,2Ddelta (t-t^{prime} ),end{array}$$

in which ({mathbb{E}}[,cdot ,]) is the expectation of the variable and D is a diffusion constant. Then, the time evolution of the defect positions is represented by the following Langevin equation:

$$dot{{boldsymbol{w}}}=-frac{1}{zeta }{nabla }_{{boldsymbol{w}}}F+{boldsymbol{R}}(t),$$

(6)

where ζ is a friction constant and ({nabla }_{{boldsymbol{w}}}={left[frac{partial }{partial {u}_{1}},frac{partial }{partial {v}_{1}},frac{partial }{partial {u}_{2}},frac{partial }{partial {v}_{2}},right]}^{T}). Note that the random force R(t) represents the randomness of the cell activity during cell proliferation and self-alignment. From Equation (6), we can derive the following Fokker–Planck equation, which describes the time evolution of the probability density function of the defect positions p(w₁, w₂):

$$frac{partial p({w}_{1},{w}_{2})}{partial t}=frac{1}{zeta }{nabla }_{{boldsymbol{w}}}cdot (({nabla }_{{boldsymbol{w}}}F)p({boldsymbol{w}}))+D{nabla }_{{boldsymbol{w}}}^{2}p({boldsymbol{w}}).$$

(7)

By setting the left-hand side of the above equation to zero, the equilibrium distribution p_eq(w) can be obtained as follows:

$$begin{array}{l}{p}_{{rm{eq}}}({w}_{1},{w}_{2})=frac{exp (-F/T)}{intintintint exp (-F/T)d{u}_{1}d{v}_{1}d{u}_{2}d{v}_{2}}\qquadqquadquad =frac{1}{Z(T)}exp left(-frac{1}{T}cdot frac{1}{2i}{oint }_{partial g^{-1}(tilde{Omega})}frac{partial {f}_{Omega }(g(z))}{partial z}overline{{f}_{Omega }(g(z))}dzright),end{array}$$

(8)

where T = 2ζD/K and

$$Z(T)=intintintint exp (-F/T)d{u}_{1}d{v}_{1}d{u}_{2}d{v}_{2}.$$

The parameter T is called the effective temperature in the literature¹³. Thus, the equilibrium distribution of the defect fluctuations can be explicitly calculated from (8) for any simply connected domain.

To calculate the equilibrium distribution (8), we should determine the parameter T, which represents the magnitude of the fluctuations caused by various biological and physical factors such as cell types and activities. We estimated T from the data shown in Fig. 1c by the maximum likelihood method as follows.

Let ({w}^{(k)}={[{w}_{1}^{(k)},{w}_{2}^{(k)}]}^{T}) be a vector composed of defect positions for the kth experimental data. We denote a set of defect positions obtained from N microwells in the experiments as ({mathcal{W}}:={{w}^{(1)},cdots ,,{w}^{(N)}}). When the trials are independent, the likelihood function (p({mathcal{W}},| ,T)) is calculated as

$$begin{array}{ll}p({mathcal{W}},| ,T),=,mathop{prod }limits_{k=1}^{N}p({w}^{(k)},| ,T)\qquadqquadquad=,mathop{prod }limits_{k=1}^{N}frac{1}{Z(T)}exp left[-frac{1}{T}F({w}^{(k)})right].end{array}$$

(9)

To get the maximum likelihood estimation of T, we should maximize the following log-likelihood with respect to T:

$$log p({mathcal{W}},| ,T)=-mathop{sum }limits_{k=1}^{N}frac{F({w}^{(k)})}{T}-Nlog Z(T).$$

(10)

By changing the variable T to the inverse temperature β: = 1/T, (log Z(beta )) can be calculated by a numerical integration of the following equation according to ref. ²⁷:

$$log Z(beta )=log Z(0)+mathop{int}nolimits_{!0}^{beta }{mathbb{E}}[-F({w}^{(k)})]dbeta ,$$

(11)

where Z(0) is the square of the area of the domain Ω and ({mathbb{E}}[-F({w}^{(k)})]) is the expectation of −F(w^(k)), which can be calculated by an MCMC method. Figure 4a shows the log-likelihood function. The log-likelihood function had a single peak, which corresponded to the maximum likelihood estimation of T.

**Fig. 4: Estimation of the magnitude of the stochastic fluctuation based on the MCMC method.**

Figure 4b shows the estimated parameters T, which did not vary even when α increased. From the estimation, the spatial distributions of the defect positions were sampled by MCMC simulations (Fig. 5). In the result, the spatial distributions of the simulated defects were localized around the horizontal axis, indicating that the obtained densities were qualitatively the same as the experimental results (Fig. 1c). This consistency supports the validity of the proposed estimation method for the parameter T.

Using the experiment results (Fig. 1c) and MCMC simulations (Fig. 5), histograms of the defect positions and simulated distributions of defects were calculated (Fig. 6). The experimental and numerical results agreed in the sense that the histograms of the absolute value of the u-component of the defect position were right-biased (Fig. 6a) and that the histograms of the angle between the defect and the u-axis were symmetric (Fig. 6b). Moreover, the distributions of (theta) narrowed as α increased. Therefore, it was verified that the proposed stochastic model was consistent with the experimental results.

**Fig. 6: Histograms of the defect positions and simulated distributions of defects.**

It is noted that the peaks of the distributions of u for the experimental data (blue histograms in Fig. 6a) was slightly shifted to the left compared to those for the MCMC simulations (red lines in Fig. 6a), which means the defects tended to be generated closer to the origin than the theoretical positions. A similar shift was also reported for C2C12 cells confined in circular geometries of radius 200 to 400 μm¹³. Since some research has reported that C2C12 cells show contractile behavior as active nematic liquid crystals¹⁷, this result suggests that we should consider the active contractile forces as well as the forces calculated from the Frank elastic energy for asymmetric geometries. However, since the proposed method can precisely calculate defect positions without active stresses, we may be able to estimate the magnitude of the active stresses from the gap of the peaks between the experimental and theoretical histograms by culturing cells in asymmetric geometries. Our future work will include solving such inverse problems of the active stresses of aligned cells in asymmetric geometries from experimental data.

Discussion

In the proposed stochastic model of the defect positions in non-circular geometries based on the Fokker–Planck equation, the parameter T estimated from the data is considered to be the degree of the stochastic fluctuations of the cell activities. Duclos et al.’s study for circular geometries¹³ estimated T to be 0.1 to 0.2, and they concluded that the parameter T does not significantly depend on the size of confinement and cell types. However, they could not investigate the geometric effects of local curvature of the confinement on the defect positions and their randomness, because their proposed model was confirmed only for circular domains. In contrast to their study, our proposed method presented in this paper provides a calculation method for estimating the parameter T for any simply connected domains. Although the estimated value of T was slightly larger in this study, it was still smaller than 1, which suggests our result has the same trend as Duclos et al’s study. Furthermore, T was almost the same for different shapes, indicating that the stochasticity of the defect distribution was not affected by the geometry. Thus, our findings demonstrated that the cell activity for self-alignment was robust against changes in confinement geometry, and thus, the behavior of the cell alignment in non-circular domains could be well predicted using the proposed theory, which is based on nematic liquid crystal theory. From the viewpoint of biology, we think these findings should contribute to a better understanding of the alignment order and behavior of biological cells and provide a mathematical method for the predictive design of cellular tissues for experiments. Specifically, since our approach can be applied to any bounded and simply connected domain, the proposed method offers a precise and quantitative analysis of stochastic fluctuations of topographical defects in various kinds of geometries, which will allow the optimization and predictive design of cell alignment in patterned cell populations for experiments.

The proposed prediction method for defect distributions can be applied to any simply connected domain since the cell alignment and its Frank elastic energy are explicitly expressed by complex potentials and conformal mappings according to our previous study²⁵. Although the present study focused only on simply connected domains, the same method should be applicable to a different topology if we can explicitly express the Frank elastic energy with respect to the defect positions and conformal mapping for the topology. At the very least, we know that we can apply the proposed method to an annulus or any doubly connected domains because our group recently obtained an explicit formula for cell alignment in doubly connected domains²⁸. We also expect the proposed method can be applied to any multiply connected domains by constructing the explicit formula for those domains, which is our future work.

This study proposed a method to quantify stochastic fluctuations of topological defects in cell populations confined to a simply connected domain. We modeled the stochastic fluctuations of topological defects in confined cell populations as random walks in external forces calculated from the Frank elastic energy. We then showed that the equilibrium distributions of the defect positions could be explicitly calculated using conformal mappings and MCMC simulations. The proposed method was verified by comparing it with the results of C2C12 cells cultured on microwells. The experimental results and numerical simulations both showed that the defect positions were localized as the geometries became constricted. The estimated temperature parameter T was small and did not vary when the geometry changed, indicating that the spatial order of the cell alignment can be predicted by the proposed stochastic model and the estimated temperature parameter T.

Methods

Markov chain Monte Carlo simulations

To obtain the spatial distributions of the defect positions from Equation (8), we used the Metropolis–Hasting method for Gibbs canonical distributions²⁹.

Let ({W}^{(k)}:={{w}_{1}^{(k)},{w}_{2}^{(k)}}) be the set of sampled defects for the kth trial. The elastic energy for the positions W^(k) is defined by F(W^(k)). To sample the k + 1th defect, we determine the set of the next candidate samples (W^{prime}) as follows. We randomly choose one defect ({w}_{m}^{(k)}) (m = 1 or 2) from the set W^(k). The next candidate position for the m − th defect is then randomly chosen from the intersection of the domain Ω and a tiny square with vertices (({rm{Re}}[{w}_{m}]pm delta ,{rm{Im}}[{w}_{m}]pm delta )). For the other defect, the next candidate position is the same as the kth trial. In this way, we can obtain the next candidate (W^{prime}) and calculate its elastic energy (F(W^{prime} )). By comparing (F(W^{prime} )) with F(W^(k)), we determine the set of the next samples W^(k+1) as follows:

1.

If (F(W^{prime} )le F({W}^{(k)})), ({W}^{(k+1)}=W^{prime}).
2.

Otherwise, we sample a random number ξ from the uniform distribution on (0, 1) and determine W^(k+1) as follows:

$${W}^{(k+1)}=left{begin{array}{ll}W^{prime} ,quadquad xi le exp (({W}^{(k)}-W^{prime} )/T)\ {W}^{(k)},quad; {rm{otherwise}}end{array}right.$$

In the MCMC simulations, the observed defect positions were normalized so that 300 μm was equal to 1. For the calculation of the elastic energy (5), ϵ was set to 0.01. When choosing the next candidate position in the MCMC simulations, δ was set to 0.1. The range of the inverse temperature β was 0 ≤ β ≤ 4, which was divided into steps of 0.1 for the numerical integration in Equation (11), and the trapezoidal rule was used. To obtain stationary distributions independent of the initial values, we take 2 × 10⁵ samples and truncate the first 10³ samples according to the Gelman–Rubin method³⁰.

Fabrication of cell culture microwells

To make polydimethylsiloxane (PDMS) microwells for cell culturing, SU-8 molds were fabricated on silicon (Si) substrates through the following steps. SU-8 dry film with a thickness of 45 μm (SU-8 3045CF DFR (WF1), Nippon Kayaku, Japan) was laminated on a Si substrate. Then, UV light was illuminated onto the SU-8 film through a photo mask. The substrate was pre-baked on a hot plate for 5 min at 100 °C. After the prebake, the substrate was immersed in SU-8 developer (Kayaku Advanced Materials, USA) to develop the resist and rinsed with the developer. Finally, the substrate was hard-baked on a hot plate for 1 h at 180 °C.

PDMS microwells were fabricated as follows. First, a base elastomer and curing agent of SILPOT184 (DuPont Toray Specialty Materials, Japan) were mixed in a weight ratio of 10:1 and degassed. The mixture was poured onto the SU-8 mold and baked in an oven for 3 h at 60 °C. After the bake, the PDMS substrate was peeled from the mold and cut with a scalpel knife.

Before culturing cells on the PDMS microwells, the microwells were coated with adhesion proteins. First, the substrate was sterilized by irradiating UV light in an incubator for 1 h. Then, the surface of the substrate was hydrophilized by vacuum plasma (YHS-R, SAKIGAKE-Semiconductor, Japan) for 5 s and coated with a 5 μg/mL fibronectin solution (063-05591, Fujifilm Wako Pure Chemical Corporation, Japan) in phosphate-buffered saline (PBS(–), 14249-95, NACALAI TESQUE, Japan) for 30–60 min. Finally, the substrate was rinsed with PBS(–) buffer.

Cell cultures

Mouse myoblast cells (C2C12) were provided by the RIKEN BioResource Center. We prepared a growth medium and a differentiation medium. The growth medium was composed of Dulbecco’s Modified Eagle Medium (DMEM, 08458-45, NACALAI TESQUE, Japan) supplemented with 10% fetal bovine serum (S-FBS-NL-015, Serana Europe GmbH, Germany) and 1% penicillin-streptomycin solution (26253-84, NACALAI TESQUE, Japan). The differentiation medium was composed of DMEM supplemented with 2% horse serum (26050070, Thermo Fisher Scientific, USA) and 1% penicillin-streptomycin solution.

The cells were cultured in an incubator at 37 °C and 5% CO₂. When passaging the cells, trypsin-EDTA solution (26253-84, NACALAI TESQUE, Japan) was used to detach the cells and prepare cell dilution with the growth medium. The cells were seeded on the PDMS substrates and cultured in the growth medium for 1 to 2 days to a confluent state. Then, the growth medium was replaced with the differentiation medium to suppress cell division, and cultured for 1 day. Finally, the cultured cells were observed with an inverted phase-contrast microscope (IX-70, Olympus, Japan) equipped with a 1.5 × magnification changer, a 4 × objective lens (UPlanF1, Olympus, Japan), and a digital CCD camera (ORCA-R2, Hamamatsu Photonics, Japan).

Image processing

The phase-contrast images were analyzed with FIJI³¹ as follows. Using OrientationJ³², the local alignment angle of the cells ϕ(u, v) was obtained. From ϕ(u, v), the following scalar order parameter S was calculated:

$$S=sqrt{{langle cos 2phi rangle }^{2}+{langle sin 2phi rangle }^{2}},$$

(12)

where 〈 ⋅ 〉 is the spatial average in a 20 pixel × 20 pixel (≈20 μm × 20 μm) area. Local minima of the order parameter S were considered as candidate points of topological defects. Each candidate point was checked to determine whether it is a defect by calculating the winding number around the candidate point.

Derivation of the holomorphic function h(z)

Using the conformal mapping g(z) from the unit disk ({mathbb{D}}) to the simply connected domain Ω, the boundary of the domain ∂Ω is parameterized as follows:

$$partial Omega :={w(theta )=g({e}^{itheta })=u(theta )+iv(theta ),| ,theta in [0,2pi )}.$$

The tangent vector at w(θ) = g(e^iθ) is ({left[frac{partial u}{partial theta },frac{partial v}{partial theta }right]}^{T}) and the tangential angle at w(θ) is expressed as

$$nu (theta )=arg left(frac{partial u}{partial theta }+ifrac{partial v}{partial theta }right)=arg left(frac{dw(theta )}{dtheta }right).$$

(13)

Furthermore, (frac{dw(theta )}{dtheta }) is calculated as follows:

$$frac{dw(theta )}{dtheta}={frac{dw(z)}{dz}frac{dz}{dtheta}Big|}_{z = {e}^{itheta}}=i{e}^{itheta}g^{prime} ({e}^{itheta}).$$

(14)

Here, we consider the alignment angle ϕ defined in [−π/2, +π/2) and identify ϕ +nπ with ϕ for any integer (nin {mathbb{Z}}), as nematic liquid crystals are nonpolar molecules. Then, we can consider ϕ in the quotient space ({mathbb{R}}/pi {mathbb{Z}}), for which (arg (i)=frac{pi }{2}=-frac{pi }{2}). Thus, we have

$$arg left(frac{dw(theta )}{dtheta }right)=left(theta -frac{pi }{2}right)+arg (g^{prime} ({e}^{itheta })).$$

(15)

The first term ((theta -frac{pi }{2})) is equal to the real part of the complex potential ({f}_{{mathbb{D}}}(z,| ,{z}_{1},{z}_{2})) when z = e^iθ according to Proposition 1 in ref. ²⁵. The second term is equal to the real part of (-ilog (g^{prime} ({e}^{itheta }))). Thus, the real part of the complex potential (1) is equal to the tangential angle ν(θ) at the boundary ∂Ω. Since the function g(z) is conformal, (g^{prime} (z)ne 0) for any (zin {mathbb{D}})³³, and thus the potential (-ilog g^{prime} (z)) is holomorphic in Ω. Accordingly, the complementary term h(z) is explicitly expressed as (h(z)=-ilog g^{prime} (z)).