Explicit calculation method for cell alignment in non-circular geometries.

The alignment of spindle-shaped cells in two-dimensional geometries induces singular points called topological defects, at which the alignment angle of the cell cannot be defined. To control defects related to biological roles such as cell apoptosis, calculation methods for predicting the defect positions are required. This study proposes an explicit calculation method for predicting cell alignment and defect positions in non-circular geometries. First, a complex potential is introduced to describe the alignment angles of cells, which is used to derive an explicit formula for cell alignment in a unit disc. Then, the derived formula for the unit disc is extended to the case for non-circular geometries using a numerical conformal mapping. Finally, the complex potential allows a calculation of the Frank elastic energy, which can be minimized with respect to the defect positions to predict their equilibrium state in the geometry. The proposed calculation method is used to demonstrate a numerical prediction of multiple defects in circular and non-circular geometries, which are consistent with previous experimental results.

Introduction

Spindle- and rod-shaped living materials exhibit active behaviours as nematic liquid crystals, known as active nematics [1,2]. Many active nematic systems have been studied, including micro- and nano-biosystems such as microtubule-kinesin compounds [3,4], bacterial colonies [5-8], cell monolayers [9-12], macroscopic tissues [13,14] and animals [15]. In these active nematic systems, the alignment of the spindle-shaped materials induces singular points at which the alignment angle cannot be defined. These singular points are called topological defects. The non-uniformity of the alignment around the defects can cause mechanical imbalances among the components, which may play key roles in the active movement of living matter. Therefore, many experimental and theoretical studies of active nematic systems have focused on the characteristics of the topological defects.

Recent studies have revealed that topological defects in spindle-shaped adhesive cells have different properties from other active nematic systems. Cells cultured in small domains tend to align along the boundaries of the domains. The cell alignment along the boundaries then forces the surrounding cells to align in the same direction at a scale of several hundred micrometres [16,17]. A recent study has shown that such cell alignment leads to the stable and reproducible generation of two topological defects at specific positions when the cells are cultured in a small disc [9]. The defect positions in the disc can be theoretically calculated using nematic liquid crystal theory [9]. It has also been shown that topological defects in certain cell populations can induce apoptosis, which results in extrusion of the cells [10]. Therefore, topological defects are expected to be a mechanical factor for cell function in the field of mechanobiology [18]. Moreover, new materials have been proposed for controlling the alignment and defects in biological nematic cells [19-21]. To systematically control biological functions by topological defects, it is necessary to develop a numerical calculation and design method for predicting and controlling topological defects in cell populations.

There are two approaches for calculating alignment angles of nematic systems in a domain where several topological defects exist. One uses a finite-element method (FEM) [22,23]. A recent study [24] calculated cell alignment and polarization by expressing a contraction force with nematic order parameters. The other method is to derive an explicit expression of the alignment angles in terms of the defect positions. Duclos et al. [9] derived a Fourier series expression for the alignment angles of spindle-shaped cells in a disc where two defects exist equidistant from the centre. They determined that there is a certain defect position which minimizes the Frank elastic energy by differentiating the energy with respect to a variable related to the defect positions. The predicted position was consistent with cell culture experiments. Therefore, explicit expressions of the cell alignment can reveal the spatial dependence of the defect positions on the total elastic energy. However, there are some limitations in the explicit expression derived in the previous study; it is based on a Fourier series expansion on a disc and cannot be applied to non-circular geometries. In addition, they considered only a special case when there exist two, symmetrically spaced defects in a disc. To predict defect positions in non-circular geometries, it is necessary to derive general expressions for the cell alignment in non-circular geometries, such as the star-shaped geometry reported in Saw et al. [10]

In this paper, we propose an explicit method for calculating cell alignment in non-circular geometries. The proposed method employs a complex-analysis-based approach to express the cell alignment, which enables (i) a simpler expression compared with a Fourier series expression, (ii) generalization of the geometries using conformal mappings, and (iii) easy computation of the elastic energy using contour integration. This paper is organized as follows. In §2.1, we review a physical model and a Fourier series solution for the alignment angles of spindle-shaped cells in a disc according to Duclos et al. [9]. Then, we derive an explicit expression for the cell alignment in a unit disc using complex potentials in §2.2. In §3, we combine the derived explicit expression with conformal mapping techniques to compute the cell alignment in non-circular geometries. In §4, we propose a numerical minimization method for the elastic energy by calculating the contour integrals and taking a derivative of the elastic energy with respect to the defect positions. In §5, the proposed method is confirmed by several numerical simulations. Finally, we compare the simulation results with previous experimental results.

In this paper, we use the following mathematical notation. and are the sets of integers and real numbers, respectively. denotes the complex plane with coordinate z = x + iy (). is the open unit disc in , and its boundary is defined as . We denote the closure of as . For a certain domain whose boundary ∂S is a Jordan curve, we denote unit tangential and normal vectors at the boundary as and , respectively. For a scalar-valued function ϕ(x, y) on the domain S, we define tangential and normal derivatives of ϕ at the boundary as and , respectively.

Explicit expression of cell alignment in a unit disc

Physical model and Fourier series solution of cell alignment in a disc

In this section, we review a physical model and a Fourier series solution of the cell alignment in a disc according to Duclos et al. [9]. We model the spindle-shaped cells as rod-like particles with negligible size, aligned along a specific axis denoted by the unit vector d = (cosϕ(x, y), sinϕ(x, y))T, where ϕ(x, y) is the alignment angle of a rod, which varies continuously with respect to the position (x, y). The vector d is called the ‘director’ [25]. We also assume that d and −d are indistinguishable; that is, is identified with ϕ(x, y). We define a restricted value of ϕ(x, y) for [ − π/2, + π/2) as follows:where is a floor function defined in . Although both ϕ(x, y) and express the same direction, the two functions are distinguished in our paper.

Let us consider a special case where the cells exist in a disc . On the boundary ∂D, the cells are forced to align along the disc. That is, the following boundary condition is assumed:where θ denotes the polar angle on the boundary.

To derive the governing equation of the angle ϕ(x, y), Duclos et al. assumed the physical condition called the ‘one constant approximation’ [25]. Then, for a subdomain S ⊂ D that does not contain singular points corresponding to topological defects of the cells, the total elastic energy (i.e. Frank elastic energy) F is expressed as follows:where K is the Frank elastic constant. By taking the variation of F with respect to ϕ(x, y), it is easily found that the angle ϕ(x, y) should satisfy the Laplace equationTherefore, ϕ(x, y) is a harmonic function.

For a precise definition of topological defects, Duclos et al. introduced another scalar function ψ(x, y) which satisfies , where = [0, 0, 1]T. It is noted that this constraint is equivalent to the Cauchy–Riemann equations and thus ψ(x, y) is the conjugate harmonic function of ϕ(x, y). Using ψ(x, y), a singular point is defined as a point that satisfies the following Poisson equation:Let be a bounded domain that contains a single singular point ; it then holds from Stokes’ theorem and equation (2.5) thatwhere d is a line element along . This equation implies the angle of a cell changes by +π around the singular point in the counter-clockwise direction. We call this singular point a topological defect at with a topological charge of +1/2 (or a +1/2 defect). More generally, we define a topological defect at with a topological charge of a half-integer q = m/2 () as a singular point which satisfies the following Poisson equation:In the same way as in equation (2.6), we can show that the angle changes by 2qπ around the defect with a topological charge of q.

In the study of Duclos et al., it was found experimentally that only two +1/2 defects existed when spindle-shaped cells were cultured in a disc D. To obtain the distribution of the angles of the cells in the disc, they assumed that the two +1/2 defects were placed at the same distance r0 from the centre of the disc and considered the following Poisson equation:whereNote that it follows from the boundary condition (2.2) and the Cauchy–Riemann relation ∂ϕ/∂s = −∂ψ/∂n that ∂ψ/∂n = −1 on the boundary ∂D. Under this boundary condition for ψ, they solved the Poisson equation (2.8) using the Fourier series of ψ. Considering the Cauchy–Riemann relations between ϕ and ψ, the following closed form of ϕ was obtained.

Explicit expression of cell alignment using complex potentials

As described in §2.1, the previous study [9] modelled the cell alignment as a nematic liquid crystal in a disc and derived the closed form of ϕ for a special case when there exist two, symmetrically placed, topological defects. In this subsection, we give a simpler and more general expression of the alignment angle ϕ(x, y) for arbitrary positions and number of defects in the unit disc. To this end, we re-formulate the problem described in §2.1 using complex potentials.

We assume that N topological defects exist at z = x + iy, where |z| < 1, for k = 1, …, N. We denote the topological charge of the kth defect by q for k = 1, …, N and the set of defect positions by . Let be the restricted value of ϕ(x, y) to [ − π/2, + π/2) calculated by equation (2.1). Considering equations (2.2), (2.4) and (2.6), the alignment angle of the cell ϕ(x, y) satisfies the following three conditions: The following proposition provides an explicit expression of ϕ(x, y) which satisfies the above three conditions.

ϕ(x, y) is a harmonic function in .

For any simply connected, bounded domain D ⊂ D which contains only the kth defect, ϕ(x, y) satisfies

For the boundary , is equal to its tangential angle. That is,

Proposition 2.1.

For a holomorphic function in given bywhereits real part, , satisfies conditions (i)–(iii).

Proof.

Condition (i) is clearly satisfied because log (z − z) and are holomorphic on and , respectively, and hence their real parts are harmonic.

For condition (ii), for the kth domain D, since is harmonic inside the domain D for l = 1, …, N, we have from Green’s theoremIn the same way as described above, when l ≠ k, because arg(z − z) is harmonic inside D. Thus, we need to consider only the term . Let D(ε) : = {z + r ei | 0 ≤ r ≤ ε, 0 ≤ θ < 2π}, where ε is a small positive number. Since arg(z − z) is harmonic on , it holds from Green’s theorem thatThen, it holds thatFrom the above discussion, we haveHence condition (ii) holds for .

For condition (iii), on the boundary , is equal to up to a constant 2Lπi, where L is an integer. Therefore, the real part of the sum of the pair of and becomesSince q is a half-integer, the restriction of equation (2.19) to [ − π/2, π/2) by (2.1) becomes qθ. Using equation (2.14), the restricted value satisfies condition (iii). ▪

Remark 2.2.

The condition in proposition 2.1 is necessary for the proof of condition (iii). Conversely, we can show is equal to 1 when satisfies conditions (i)–(iii) as follows. The total change of along the boundary should be equal to 2π from condition (iii). However, because is harmonic in , the line integral is equal to from condition (ii) and Green’s theorem. Therefore, should be equal to 1.

Our expression for in proposition 2.1 is equivalent to equation (2.10) derived in the previous study [9], when N = 2, q1 = q2 = 1/2, and the defect positions are given by equation (2.9). Moreover, the imaginary part of , denoted by , for that case satisfies the Poisson equation (2.8) because −(1/2π)log|z| is a Green function of the two-dimensional Laplace operator. Therefore, our result is a general expression for the angle ϕ(x, y) as compared with equation (2.10) without any assumptions on the defect number and positions in the unit disc.

It is noted that our expression can be obtained by the method of images in electromagnetics. in represents a complex potential by the pole at a mirror image of the original kth defect with respect to the unit circle . In other words, is the potential of a defect with the same topological charge q located at (the reflection of the point z about the unit circle). Therefore, the derived potential comprises the potentials of the original defects, log (z − z), and those of their mirror images, . Figure 1 shows an example of where N = 2, q1 = q2 = +1/2, z1 = 0.75ei, and z2 = 0.85e−i. The alignment angle is exactly equal to the tangential angle of the boundary of the unit disc (figure 1a) and increases linearly in [0, π) and [π, 2π) along the unit circle (figure 1b).

Figure 1

Calculation of cell alignment when two +1/2 defects exist in the unit disc. (a) Colour map of . The blue and green circles indicate the original defects and their mirror images about the unit circle, respectively. (b) Boundary values of and along the unit circle. We computed from by equation (2.1).

Explicit calculation of cell alignment in non-circular geometries using conformal mappings

In this section, we extend the explicit expression (2.13) to the generalized case for a simply connected and bounded domain in which is bounded by a Jordan curve as shown in figure 2. We denote points of the domain and the open unit disc by w = u + iv and z = x + iy, respectively, which means that and are defined on the w-plane and the z-plane, respectively. We assume that there are N topological defects at w = u + iv with a topological charge of q for k = 1, …, N. Considering proposition 2.1 and remark 2.2, we assume the sum of the topological charge to be equal to 1. We denote the set of defect positions in Ω by . According to the Riemann mapping theorem, there always exists a one-to-one conformal mapping g that maps the open unit disc to [26]. We denote the set of topological defects mapped from . Since the boundaries and are Jordan curves, it follows from Carathéodry’s theorem that the conformal mapping g extends continuously and one-to-one to the closed unit disc [27]. Therefore, the boundary of the domain is parametrically defined as the image of using g

Figure 2

Schematic of geometries and defect positions for the calculation of cell alignment in a non-circular geometry .

The problem here is to calculate the alignment angle ϕ(u, v) under the condition that the cells align along the boundary . Using the same discussion as in §2.1, we can show that the alignment angle ϕ(u, v) should satisfy the Laplace equation in a subdomain of which does not contain any topological defects. In addition, a topological defect in is defined as a singular point of the conjugate harmonic function ψ(u, v) of ϕ(u, v) which satisfies equation (2.7). Therefore, ϕ(u, v) should satisfy the following three conditions:

ϕ(u, v) is a harmonic function on .

For any simply connected, bounded domain D which surrounds only the kth defect, ϕ(u, v) satisfies

On the boundary , the restricted value is equal to its tangential angle ν(θ):

Remark 3.1.

ν(θ) is expressed in terms of the conformal mapping g as follows. Since the boundary is a parametric curve (3.1), the arc-length of between g(ei0) and g(ei) is given bywhich gives the following equation:According to the theory of elementary differential geometry [28], the curvature of the curve κ is expressed as κ = dν/ds. Therefore, ν(θ) is expressed in terms of the curvature asAccording to complex analysis [29], the curvature is calculated from the conformal mapping g(z) as

Let us derive an explicit expression of ϕ(u, v) which satisfies conditions (I)–(III). To this end, for ν(θ), we defineThen, the Poisson–Schwartz integral formula giveswhich is holomorphic in and whose real part on coincides with μ(θ). Now, by using h(z) and the conformal mappings g(z) and g−1(w), an explicit expression of ϕ(u, v) which satisfies conditions (I)–(III) is given as follows.

Theorem 3.2.

Let be the following holomorphic function on :where is defined by equation (2.13) for z = g−1(w) and . Then, its real part, , satisfies conditions (I)–(III).

Proof.

We can consider that is a composite function of the conformal mapping w = g(z) and can be viewed as a function defined on the unit disc . In what follows, we change variables from w to z and consider conditions (I)–(III) on the z-plane. To this end, we denote the real part by

We first show that satisfies condition (I). Since g(z) and g−1(w) are holomorphic,According to proposition 2.1, is harmonic as a function of x and y on . In addition, Re[h(z)] is also harmonic on , from the Poisson–Schwartz formula. Therefore, the right-hand side of equation (3.12) is zero on and hence the sum of these real parts, , is harmonic on , which means condition (I) holds.

For condition (II), it holds thatSince g−1(w) is a one-to-one mapping and D encloses only the kth defect located at w, the inverse g−1(D) in encloses only the kth defect located at z = g−1(w). Moreover, satisfies condition (ii) in §2.2. Therefore, the first contour integral in the right-hand side of equation (3.13) is equal to 2πq. For the second integral in the right-hand side of equation (3.13), since h(z) is holomorphic in , its real part Re[h(z)] is harmonic in g−1(D). Then, it follows from Green’s theorem thatThus, satisfies condition (II).

For condition (III), we substitute w(θ) = g(ei) into equation (3.10)According to proposition 2.1, is equal to θ − π/2. In addition, Re[h(ei)] is equal to ν(θ) − θ + π/2 from the definition of h(z). Therefore, the restricted value is equal to ν(θ) on the boundary . ▪

Based on theorem 3.2, we propose a method to calculate the alignment angle in a non-circular domain using conformal mappings as follows. We illustrate an example when g(z) = z + 0.2z3 in figure 3. First, the boundary condition for is interpreted as the boundary condition for on . Then, the boundary condition is decomposed to those for the real parts of and h(z) (figure 3a). Since μ(θ) (i.e. the boundary condition for Re[h(z)]) represents the difference of the two boundary conditions (2.12) and (3.3), h(z) is interpreted as an additional term for satisfying the boundary condition (3.3), which is determined only by the conformal mapping g. From the two boundary conditions, and h(z) are determined from equations (2.13) and (3.9), respectively. By adding these two functions and using the conformal mapping z = g−1(w), we can calculate the alignment angle at arbitrary points in (figure 3b).

Figure 3

Schematic view of the calculated orientational angle in a non-circular region using conformal mapping g(z) = z + 0.2z3. (a) Boundary condition for the real part of and its decomposition to the boundary conditions for the real parts of and h(g−1(w(θ)). (b) Real parts of the complex functions , , and h(g−1(w)).

Prediction of defect positions by numerical minimization of Frank elastic energy

The alignment angle expressions in §§2.2 and 3 are functions of the defect positions z or w, where those positions can be arbitrary in the domain. However, experimental studies of cell culturing have shown that defects tend to be generated on certain points which minimize the total elastic energy because cells could actively migrate [9,10]. Therefore, it is possible to predict defect positions by minimizing the elastic energy with respect to them. This section provides a systematic method for minimizing the elastic energy using formulae of complex analysis.

Since the complex potential (3.10) has singular points at defects, we should exclude the neighbourhoods of the defects in computing the elastic energy so that it does not diverge to infinity. For that purpose, we define a small disc centred at the kth defect as D(ε) = {w + r e−i| 0 ≤ r ≤ ε, 0 ≤ θ ≤ 2π} (0 < ε ≪ 1). Then, the elastic energy F is calculated aswhere . We take K/2 = 1 because the constant K is not concerned with the minimization of the energy with respect to defect positions. By using the Cauchy–Riemann equations for the complex potential (3.10), the integrand in (4.1) can be expressed by complex differentiation as follows:where and . Hence, applying the complex form of Green’s theorem , the elastic energy F can be converted to a contour integral form:Here, equation (4.3) is rewritten as a contour integral on the z-plane aswhere and . Using equation (3.10), equation (4.4) is rewritten aswhere z = g−1(w). In this way, the elastic energy can be explicitly described as a function of the defect positions z. In what follows, we consider the minimization of the energy F in with respect to z instead of the minimization with respect to w. Once the minimizer of F is obtained, w is obtained as w = g(z).

To minimize F, we differentiate F with respect to , givingSince is considered the gradient of F [30], the total energy F can be numerically minimized by the steepest descent method as follows:where α is a small positive constant.

To evaluate the contour integral of the elastic energy on , we set the integral path as shown in figure 4. We define the half-line in the z-plane from the centre of the unit disc through the defect position z as for k = 1, …, N. Let c be an intersection point in the w-plane between and and d be the intersection point in the w-plane between g(Γ) and D which satisfies |g−1(d)| > |z|. Then, we define a straight path in the z-plane from g−1(c) to g−1(d) as and its inverse given by . We denote by the contour from g−1(c) to g−1(c) in . c is identified with c1. The contour integral was evaluated along , g−1(∂D) in the clockwise direction, and in the counter-clockwise direction, for k = 1, 2, …, N.

Figure 4

Schematic diagram of integral path on the original space and the unit disc .

Numerical experiments

Computational implementation

In all the calculations in this paper, we used Matlab 2021a. To calculate the complex logarithm log (z − a), we re-defined the log function asBy this re-definition, we chose a branch cut of log (z − a) as {a + texp (i arg(a)) | t > 0}.

For numerical calculation of the cell alignment and integration of the elastic energy, we divided the boundary of the unit disc into M = 4000 equal arcs by M points ζ = ei2 (m = 0, 1, …, M − 1). We used the trapezoidal rule in calculating the contour integral in equations (4.5) and (4.6).

To compute the value of the analytic function h(z) on in theorem 3.2 from its known real part on the boundary, we used a Fourier expansion series as follows. Since μ(θ) is a real and periodic function of θ ∈ [0, 2π), μ(θ) is expressed as the following Fourier series expansion:where a and b are the Fourier series coefficients for the cosine and sine terms, respectively. By truncating the summation in equation (5.2) at the M/2th term, we have an approximated value of μ(θ) at θ = 2πm/M, given bywhere M′ = M/2. The holomorphic function h(z) on can be approximately expanded as a power series of for 0 ≤ r < 1 by truncating the series at the M′th termwhere c = d + ie (). Since the limit of the real part of equation (5.4) should continuously approach μ(θ) when r goes to 1, we have a = d and b = −e. Therefore, the value of for 0 ≤ r ≤ 1 can be calculated from the Fourier series coefficients of μ(θ), given byWe calculated the coefficients a and b from the fast Fourier transform (FFT) of μ(θ). For faster calculation of the obtained h(z), h(z) was approximated by an AAA algorithm [31] programmed in Chebfun [32]. We also used conformal in the Chebfun package to get the inverse of the conformal mapping.

For the contour integrals in (4.5) and (4.6), the radius of the disc D was set to ε = 10−2. Each integral path (, , and ∂D) was divided into 1000 equal sections. To calculate the line integral on the paths and , we multiplied by exp (iε2) and exp ( − iε2) (ε2 = 10−10) for the points on and , because the paths and are on the branch cut of the complex logarithmic function defined in (5.1).

When the energy was minimized by the steepest descent method shown in (4.7), the step size α was set to 10−3 throughout the simulations. If the distance between two defects became less than 2ε, the update for the two defects was stopped and the two defects were removed from the calculations to prevent large oscillations of the defect positions. It is noted that this event can be viewed as annihilation of a defect pair.

Experimental settings

To verify the prediction of the defects by minimizing the elastic energy, we considered the following domains: (A) the unit disc and (B) non-circular domains (p = 2, 3, 4), whose boundaries were given bywhere β was set to 0.3 unless otherwise noted.

For these domains, we performed numerical demonstrations for the following cases.

Case A: Unit disc. First, we considered the simplest case of two +1/2 defects at z1 = x and z2 = −x (0 < x < 1) in the unit disc to verify the consistency with the previous study [9] (Case A-1). The minimizer of the elastic energy with respect to x was determined by calculating the elastic energy and its derivative for various parameters ranging from 0.020 to 0.980 in steps of 0.001. Next, to show the validity of the steepest descent method (4.7), the elastic energy was minimized by the steepest descent method (4.7) (Case A-2). At the initial state, two +1/2 defects are located at z1 = 0.9ei and z2 = 0.5 in . In addition, we considered minimization by the steepest descent method for the case of three +1/2 defects at z1, z2 and z3 and one −1/2 defect at z4 (Case A-3). In this case, we put z1 = 0.4ei, z2 = 0.5ei3, z3 = 0.6ei5, and z4 = 0.7ei7 initially.

Case B: Non-circular domains We predicted the defect positions in non-circular domains (p = 2, 3, 4) by the steepest descent method (4.7) (Case B-1). The initial defect positions were set as follows. For , two +1/2 defects were put at z1 = 0.6ei and z2 = −z1. For , three +1/2 defects were put at z1 = 0.8ei, z2 = 0.8ei7, z3 = 0.8ei13 and one −1/2 defect was put at z4 = 0.01e−i/100. For , four +1/2 defects were put at z1 = 0.8ei, z2 = 0.8ei5, z3 = 0.8ei9, z4 = 0.8ei13, and one −1 defect was put at z5 = 0.01e−i/100. To demonstrate that there were multiple local minima of the elastic energy for the domain , minimization of the elastic energy was performed by the steepest descent method (4.7) for two other initial conditions: z1 = g−1(0.3), z2 = 0.95ei2, , z4 = 0 (Case B-2) and z1 = g−1(0.4), z2 = 0.95ei2, , z4 = 0 (Case B-3). To explain the multiple local minima, we calculated the energy landscapes for specific defect numbers and positions as follows. First, we considered two +1/2 defects in and parametrized the defect positions as z1 = r1ei2 and z2 = r1ei4 (r1 > 0) (Case B-4). The minimizer of the elastic energy was obtained with respect to r1 from the calculated energy landscape. Then, we considered three +1/2 defects and one −1/2 defect in (Case B-5). In this case, two +1/2 defects were fixed at the positions obtained in Case B-4. The other +1/2 defect and one −1/2 defect were located at z3 = r2 (r2 > 0) and z4 = 0, respectively. The elastic energy landscape was calculated by changing the value of r2. For Cases B-4 and B-5, we set the values of β in equation (5.6) to 0, 0.15 and 0.30.

Results and discussion

Figure 5a,b shows the parametrically calculated elastic energy and its derivative in Case A-1. As shown in figure 5b, the elastic energy was minimized when x was 0.669, which is nearly equal to the theoretical value 5−1/4 reported in Duclos et al. [9]. We also confirmed that the theoretical value 5−1/4 could be re-derived using proposition 2.1 (see electronic supplementary material). A similar result was obtained in Case A-2, where the steepest descent algorithm was used for the minimization of the elastic energy (figure 5c–e). The elastic energy monotonically decreased as the algorithm proceeded, which resulted in convergence of the defect positions |z1| and |z2| to around 0.669. In addition, we confirmed that the defects were always collinearly placed in the final state, with the positions depending on the initial point of the defects located asymmetrically.

Figure 5

Predicted defect positions which minimize the elastic energy for two +1/2 defects in a unit disc . (a,b) Numerical calculation of the elastic energy and its derivative in Case A-1. (a) Schematic of the defect positions in for Case A-1. (b) Calculated elastic energy and its derivative with respect to x. (c–e) Minimization of the elastic energy by the steepest descent algorithm in Case A-2. (c) Trajectories of the defects through the iteration of the algorithm. The circle and square markers indicate initial and final defect positions, respectively. (d,e) Temporal evolution of (d) the absolute number of defect positions and (e) the elastic energy.

We numerically confirmed that the final defect positions depended on the initial positions because the minimization problem for the disc domain has infinitely many solutions due to the circular symmetry. Even so, the proposed method can obtain an optimized defect position because the gradient of the elastic energy vanishes once the defect positions converge to one of the minimizers.

Figure 6 shows the iterative minimization of the elastic energy in Case A-3 where there were four defects in the unit disc . In this case, a pair of +1/2 and −1/2 defects approached each other as the algorithm proceeded (figure 6a–c). After the attracted defects were removed, the remaining two +1/2 defects moved to the same points where the energy was a minimum as in figure 5. Throughout the algorithm, the energy monotonically decreased (figure 6d), which means that the minimization algorithm certainly worked even when there were initially four defects in .

Figure 6

Iterative minimization of the elastic energy in Case A-3. (a) Spatial distribution of the orientation angle during the minimization algorithm. The bluish and reddish circles indicate +1/2 and −1/2 defects, respectively. At steps 250 and 2000, the red and blue circles almost overlap because the two defects approached each other between step 200 and step 250. (b) Trajectories of the defects during the algorithm. The circle and square markers indicate the initial and final defect positions, respectively. (c,d) Temporal evolution of (c) the absolute number of the four defects and (d) the elastic energy.

The iterative minimization of the elastic energy also worked for the non-circular geometry case (Case B-1), as shown in figure 7. In all geometries, the elastic energy certainly decreased throughout the algorithm iterations, and the defects were finally located around the tips of lobes. It should be noted that the results for p = 4 were similar to one reported in the cells in a star-shaped geometry, where +1/2 defects tended to be generated around the tips of the star [10].

Figure 7

Prediction of defect positions in the various non-circular shapes . (Top) p = 2, (middle) p = 3, (bottom) p = 4. (a,b) Trajectories of defects in (a) the original shapes and (b) the unit disc mapped from the original domain. (c,d) Distribution of orientation angles at (c) the initial and (d) final states. (e) Time evolution of the elastic energy. The circle and square markers indicate the initial and final positions of defects, respectively. The bluish, red and yellow circles indicate +1/2, −1/2 and −1 defects, respectively.

The equilibrium points after the iteration depended on the initial defect positions since there can be several local minima for the elastic energy function. For example, when +1/2 and −1/2 defects were closely positioned in the domain (Case B-2), the two defects were attracted during the minimization (figure 8a,b). On the other hand, when the initial distance between the two defects was longer (Case B-3), the +1/2 defect moved toward the boundary of the domain. After the iteration, there remained four defects even after the minimization (figure 8c,d). Therefore, there are at least two equilibrium states in this shape. To further investigate the dependence of the initial defect positions on the final states for these cases, we calculated the elastic energy landscapes for Cases (B-2) and (B-3). For two defects as shown in figure 8e,f (Case B-4), there was only one minimal point of the elastic energy (figure 8g). When one +1/2 and one −1/2 point were added as shown in figure 8h,i (Case B-5), there were one local maxima and two local minimal points at r1 = 0 and r1 > 0 for β = 0.15, 0.30 (figure 8j). These analyses support there being two equilibrium states shown in figure 8a,d due to the existence of multiple local minima and that the proposed method will predict which state is the final state from the initial defect positions.

Figure 8

Demonstration of multiple equilibrium states due to differences in the initial defect positions. (a,b) Attracting (Case B-2) and (c,d) repelling (Case B-3) of +1/2 (bluish lines) and −1/2 (reddish lines) defects. (a,c) Trajectories of defect positions during the minimization. The circle and square markers indicate the initial and final positions of the defects, respectively. (b,d) Time evolution of the elastic energy during the minimization. (e–j) Calculation of elastic energy landscapes for (e–g) two defects (Case B-4) and (h–j) four defects (Case B-5). (e,h) Settings of the defect positions for (e) Case B-4 and (h) Case B-5. (f, i) Distribution of orientational angles corresponding to (e,h). (g,j) Relationship between the defect positions and the elastic energy. In the calculation of (j), the two +1/2 defects located in the left half-plane were determined from the local minima calculated in (g).

The proposed method allowed the calculation of cell alignment in non-circular shapes as shown in figures 3 and 7. Although we assumed a smooth boundary in the calculation of the cell alignment in §3, it is expected that the proposed method can also be applied to piecewise smooth shapes such as polygons by using Schwarz–Christoffel mappings. However, since such non-smooth shapes have singularities at corners, the elastic energy would diverge to infinity due to the singularities in the conformal mappings. In future work, we will examine solutions for such problems and investigate the effects of the corner singularities on the prediction of the defect positions.

The proposed method also enabled direct differentiation of the elastic energy with respect to the defect positions and the minimization of the elastic energy. This calculation was possible because the constructed analytic function h(z) does not depend on the defect positions. Moreover, the integration of the elastic energy is in contour integral forms, which makes it easier to calculate the integral around the singular point z compared with surface integral forms. These advantages of the proposed method contribute to the calculation of the defect positions in non-circular domains as shown in figures 6 and 7. Especially, the result for the four-leaves domain (figure 7) is consistent with the previous result which can control positions of +1/2 defects [10]. Therefore, the proposed method could be applied to predicting the defect positions and design of geometries for controlling defect positions. Our numerical results predicted the exact minimal point for various geometries, but the results have not been experimentally confirmed. Therefore, future work will include cell culturing experiments to verify the proposed calculation method.

Conclusion

This study proposed a systematic method for calculating the alignment angle and predicting defect positions which minimize the elastic energy in cell populations at a confluent state. Using several complex analysis techniques, we derived an explicit expression of the cell alignment angle in a unit disc in terms of the defect positions. The derived expression was extended to an angle calculation for non-circular geometries by using conformal mapping. By taking advantage of the explicit expression of the alignment angle in terms of the defect positions, we proposed a method for calculating the elastic energy and its derivative with respect to defect positions and minimizing the elastic energy by a steepest descent method. The validity of the methods was confirmed by numerical simulations. The results for the unit disc were consistent with the previous study, [9] while qualitative agreement with the previous study [10] was also confirmed for non-circular geometries with four lobes. Therefore, the proposed method will be a powerful tool for designing geometries to control defect positions.