Structure--property relationships of cell clusters in biotissues: 2D analysis

TL;DR

This study investigates the relationship between cell cluster structures and biotissue properties using a detailed cell adhesion model, deriving analytical solutions and simulating complex shapes to connect microstructure with macro-mechanical behavior.

Contribution

It introduces a comprehensive 2D cell adhesion model to analyze biotissue structures, providing analytical and numerical insights into their mechanical and morphological properties.

Findings

Derived equilibrium shape equations for symmetrical structures.

Obtained analytical expressions for Young's modulus, bulk modulus, and failure strength.

Simulated stable multicellular structures and observed symmetry-breaking phenomena.

Abstract

To insight the relationships between the self-organizing structures of cells, such as the cell clusters, and the properties of biotissues is helpful in revealing the function and designing biomaterial. Traditional random foam model neglects several important details of the frameworks of cell clusters, in this study we use a more complete model, cell adhesion model, to investigate the mechanical and morphological properties of the two-dimensional (2D) dry foams composed by cells. Supposing these structures are formed due to adhesion between cells, the equilibrium formations result from the minimum of the free energy. The equilibrium shape equations for high symmetrical structures without the volume constraint are derived, and the analytical results of the corresponding mechanical parameters, such as the Young's modulus, bulk modulus and failure strength, are obtained. Numerical simulation method is applied to study the complex shapes with the volume constraint and several stable multicellular structures are obtained. Symmetry-breaking due to the volume change is founded and typical periodic shapes and the corresponding phase transformations are explored. Our study provides a potential method to connect the microstructure with the macro-mechanical parameters of biotissues. The results also are helpful to understand the physical mechanism of how the structures of biotissues are formed.

Equations

$$E_i=\frac{1}{2}{\kappa }_i\oint {\Lambda }_i^{2}d{s}_i-\frac{1}{2}\sum _{j=1}^m{\omega }_{ij}{B}_{ij}.$$

$${\Omega }_i=E_i+{\gamma }_i\oint d{s}_i+\Delta p_i\int \int d{\sigma }_i,$$

$${\kappa }_i{\stackrel{...}{\varphi }}_i+\frac{1}{2}{\dot{\varphi }}_i^3-{\gamma }_i{\dot{\varphi }}_i+\Delta p_i=0.$$

$${\kappa }_i{\dot{\varphi }}_i^2+{\kappa }_j{\dot{\varphi }}_j^2-({\kappa }_i+{\kappa }_j){\dot{\varphi }}_{ij}^2=2\omega ,$$

$$\displaystyle\frac{d}{ds}\Big{[}\kappa_{i}\dot{\phi_{i}}+\kappa_{j}\dot{\phi_{j}}-(\kappa_{i}+\kappa_{j})\dot{\phi}_{ij}\Big{]}=0,$$

$${\kappa }_i{\dot{\varphi }}_i+{\kappa }_j{\dot{\varphi }}_j-({\kappa }_i+{\kappa }_j){\dot{\varphi }}_{ij}=0.$$

$${\dot{x}}_i=\cos {\varphi }_i,{\dot{y}}_i=\sin {\varphi }_i.$$

$$\displaystyle\Omega_{i}=\oint\Big{[}\frac{1}{2}\kappa_{i}\dot{\phi}_{i}^{2}+\Delta p_{i}x_{i}\sin\phi_{i}+\gamma_{i}\Big{]}ds_{i}-\frac{1}{2}\sum_{j=1}^{m}\omega_{ij}B_{ij}.$$

$${\Pi }_i({\varphi }_i,{\dot{\varphi }}_i,x_i,{\dot{x}}_i,{\lambda }_i)=\frac{1}{2}{\kappa }_i{\dot{\varphi }}_i^2+\Delta p_i{x}_i\sin {\varphi }_i+{\gamma }_i$$

$$+{\lambda }_i({\dot{x}}_i-\cos {\varphi }_i)-\frac{1}{2}\sum _{j=1}^m{\omega }_{ij},$$

$${\kappa }_i{\ddot{\varphi }}_i-\Delta p_i{x}_i\cos {\varphi }_i-{\lambda }_i\sin {\varphi }_i=0,$$

$$\Delta p_i\sin {\varphi }_i-{\dot{\lambda }}_i=0.$$

$$\frac{1}{2}{\kappa }_i{\dot{\varphi }}_i^2-\Delta p_i{x}_i\sin {\varphi }_i+{\lambda }_i\cos {\varphi }_i-{\gamma }_i=0.$$

$${\Pi }_1=\frac{1}{2}{\kappa }_1{\dot{\varphi }}_1^2+\Delta p_1{x}_1\sin {\varphi }_1+{\lambda }_1({\dot{x}}_1-\cos {\varphi }_1).$$

$${\Pi }_3=\frac{1}{2}{\kappa }_3{\dot{\varphi }}_3^2+\Delta p_3{x}_3\sin {\varphi }_3+{\lambda }_3({\dot{x}}_3-\cos {\varphi }_3)-\omega ,$$

$${\kappa }_1{\ddot{\varphi }}_1-\Delta p_1{x}_1\cos {\varphi }_1-{\lambda }_1\sin {\varphi }_1=0,$$

$$\Delta p_1\sin {\varphi }_1-{\dot{\lambda }}_1=0.$$

$${\kappa }_3{\ddot{\varphi }}_3-\Delta p_3{x}_3\cos {\varphi }_3-{\lambda }_3\sin {\varphi }_3=0,$$

$$\Delta p_3\sin {\varphi }_3-{\dot{\lambda }}_3=0.$$

$${\kappa }_1{\dot{\varphi }}_1^2+{\kappa }_2{\dot{\varphi }}_2^2-{\kappa }_3{\dot{\varphi }}_3^2=2\omega ,$$

$${\kappa }_1{\ddot{\varphi }}_1+{\kappa }_2{\ddot{\varphi }}_2-{\kappa }_3{\ddot{\varphi }}_3=0,$$

$${\kappa }_1{\dot{\varphi }}_1+{\kappa }_2{\dot{\varphi }}_2-{\kappa }_3{\dot{\varphi }}_3=0.$$

$${\dot{\varphi }}_1=-{\dot{\varphi }}_2,{\ddot{\varphi }}_1=-{\ddot{\varphi }}_2,{\dot{\varphi }}_3={\ddot{\varphi }}_3=0.$$

$$\omega =\kappa {\dot{\varphi }}_1^2,(s=s_0).$$

$$\ddot{\varphi }=\lambda \sin \varphi ,(0<s<s_0).$$

$${\dot{\varphi }}^2=2\lambda (C-\cos \varphi ),$$

$$ds=\frac{d\varphi }{\sqrt{2\lambda (C-\cos \varphi )}}.$$

$$x_0$$

$$y_0$$

$$s_0+\frac{B}{2}={\int }_0^{{\varphi }_0}\frac{d\varphi }{\sqrt{2\lambda (C-\cos \varphi )}}+\frac{B}{2}={\varphi }_0.$$

$$\sqrt{\frac{2}{\lambda (C-1)}}\times \text{FN}[\frac{{\varphi }_0}{2},\frac{2}{1-C}]+\frac{B}{2}={\varphi }_0.$$

$$\lambda =\frac{\omega }{2(C-\cos {\varphi }_0)}.$$

$$B=2{\varphi }_0-4\sqrt{\frac{C-\cos {\varphi }_0}{\omega (C-1)}}\times \text{FN}[\frac{{\varphi }_0}{2},\frac{2}{1-C}].$$

$$R$$

$$H$$

Full text

I Introduction

Soft stacking system which is composed by cells Skalak1981 , particles Seth as well as other soft units Buchcic ; Petroff ; Pennybacker often presents particular arrangement. Comparing the soft stacking system with the atom stacking system in solid state physics, the former has more complex structures. A reason for inducing this complexity is that, unlike the rigid spherical atoms, the basic units in soft stacking system have strong adaptability to change their shapes under different circumstances. For example, cells can change their shapes to explore their environment Taloni . Another reason is that, the composing unites in a soft stacking system possibly are multiplicate. Which implies that, when plenty of cells with different size (volume and surface area) are adhered together, the system will be very intricate. Due to the complexity, to insight the relationship between the configuration of soft stacking systems and the corresponding bio-functions still has challenge.

Adhesion which plays an important rule in forming the cell or vesicle self-organizing structures mainly derives from the interface interaction Engler ; Schwarz , protein interactionZhangX as well as charge interaction Pantazatos . At molecular level to simulate the three kinds of interactions in a system is constrained by the complexity and long time expending. So, continuous model provides a valuable option to reveal how the cell stacking system are formed due to adhesion. When the membrane proteins are ignored, cell membrane which is composed mostly by two layers of phospholipids can be taken as liquid bilayer vesicles as well as 2D surfaces. The equilibrium shapes for free vesicles are governed by the Helfrich-Canham bending energy theory Helfrich ; Canham ; Ouyang1 . Seifert and Lipowsky Seifert1 ; Lipowsky ; Seifert2 extended this theory to deal with the vesicle adhering system. They derived the equilibrium shape equations and boundary conditions by minimizing the free energy and also conveyed that equilibrium shapes result from the competition between the elastic energy and the adhesion energy. In past two decades, their theory has been generally used to investigate vesicle adhesion configurations and great progress has been achieved. It explained the adhesion shape composed by several red blood cells Ziherl . Deserno et al. Deserno developed a general geometrical framework to derive the equilibrium shape equations and boundary conditions, which makes it possible to reveal the complex structures composed by a large number of cells.

Cells in plant tissues often present regular frameworks Villalon ; Chen . Some example are shown in Fig. 1. We can see the length of epidermal cell can reach about 1 mm, but the width is no more than 40 $\rm{\mu}m$. Actually, the length of some plant cells can attain dozens of centimeters length. For example, the length of ramie cell can reach 145 mm, which is 3000 times of its width Pandey . So to study the equilibrium shapes of these tubular formations can be reduced to investigate the shapes of their cross sections in 2D case. Also, these structures be previously named 2D dry foams and the properties have been wildly investigated Weaire1 ; Weaire2 ; Schliecker . However, in these works several important details of the clusters of cells are neglected. In the solid foams, the walls of all cavities form a network and the joint points between them are rigid connections Schliecker ; Roberts ; Kraynik1 ; Vecchio , that is different to the corresponding points in Fig. 1(f). For example, the points $a$ to $e$ are the initial contact points between neighbor cells. Nearby each point, two edges of two neighbor cells have the asymptotic behavior to get close smoothly. After that, two edges will contact together until the next initial contact point. When the structure is under the external load, these initial contact points can shift along the edges and the contact length will change meanwhile. Further, in most cases, we can see in Fig. 1(c) that three neighbor cells fence a triangular space, such as the area $B$. Also, four neighbor cells can fence a quadrangular space, such as the area $A$. These spaces do not belong to any cells. Each cell in bioissues has its own edges and volume. In the contact region, the edge is made up of two layers of individual cell edges. We can also pick up a cell from biotissues without breaking any other cells. These phenomena indicate the biotissues are formed by plentiful cells stacked together. But in most of solid foam models, the edges of cells are taken as an integrated network and the unique character of each cell is neglected.

Based on the above analysis, it needs a more complete model to study the frameworks of cells. In this model, the biotissue should be built up by cells and its macroscopical characters should rely on the mechanical properties of every single cell. We think the cell adhesion theory is a more suitable option. There are at lest two reasons. First, for single cell, we have known that the shape is determined by the minimization of the Helfrich-Canham bending energy. Particularly, although the plant cells have thick cell walls, we can take them as thin elastic shells and this theory also is suitable (neglect the in-plane strain) Tu . Seconde, for plentiful cells, the cell adhesion theory not only considers the independence of singer cells but also contains the adhesion action between cells. In this work, we use the cell adhesion theory to study the mechanical and morphological properties of the cell clusters in 2D case. In the Sect. II, the basic theory is shown and the numerical simulation method as well as the software are introduced. In the Sect. III, some analytical results for periodic structures without volume constraint are derived. In the Sect. IV, we show some numerical results obtained by simulation. Finally, these results are recapped in a short discussion in Sec. V.

II Theoretical methods

II.1 The basic model

In 2D case, let us consider a system which is composed by $n$ cells adhered together. Supposing that the $i$th cell has $m$ neighbors adhered to it, its free energy is Seifert1 ; Lipowsky ; Seifert2 ; Deserno

[TABLE]

Here, $\kappa_{i}$ is the bending rigidity, $\Lambda_{i}$ is the curvature, $ds_{i}$ is the element of the arc length of the cell, $\omega_{ij}$ and $B_{ij}$ are the adhesion potential (work of adhesion) and adhesion length between the $i$th cell and its $j$th neighbor, respectively. Considering the length and volume constraints, the corresponding energy functional is

[TABLE]

where $\gamma_{i}$ is the line tension coefficient, $\Delta p_{i}$ is the osmotic pressure difference between the inside and outside of the cell, $d\sigma_{i}$ is the element of the area enveloped by the cell. For the whole system, the total energy is $E_{t}=\sum_{i=1}^{n}E_{i}$ and the energy functional is $\Omega_{t}=\sum_{i=1}^{n}\Omega_{i}$. The equilibrium shapes are determined by the minimization of $\Omega_{t}$.

Supposing the whole system is in the $x-y$ plane. For the $i$th cell, let $\phi_{i}$ be the angles between the tangent of each arc and the $x$ axis, and define the clockwise direction as the positive direction of $\phi_{i}$, there is $\Lambda_{i}=d\phi_{i}/ds_{i}=\dot{\phi}_{i}$. The equilibrium shape for whole system needs the first variation $\delta\Omega_{t}=0$. Due to each cell has independent length $L=\oint ds_{i}$ and volume $V_{i}=\int\int d\sigma_{i}$, $\delta\Omega_{t}=0$ yields $\delta\Omega_{i}=0$ which gives the general shape equation Ouyang1 ; ZhangS1 ; ZhangS2

[TABLE]

Further, at the initial adhesion points between $i$th cell and $j$th cell, $\delta\Omega_{t}=0$ yields the boundary conditions Deserno

[TABLE]

Considering $\phi_{i}$ and $\phi_{j}$ turn clockwise in the $i$th cell and $j$th cell, respectively. In the adhesion region we define $\phi_{ij}$ belongs to the $i$th cell ($\phi_{ij}$ turns clockwise in the $i$th cell) and $\phi_{ji}$ belongs to the $j$th cell ($\phi_{ji}$ turns clockwise in the $j$th cell) and we have $\dot{\phi}_{ij}=-\dot{\phi}_{ji}$.

There is another expression for the Eq. 3 by using the Lagrange equation. Considering $\dot{\phi}_{i}=d\phi_{i}/ds_{i}$, there are

[TABLE]

Making use of $dV_{i}=x_{i}\sin\phi_{i}ds_{i}$, Eq. 2 is reduced to

[TABLE]

Consequently, the Lagrange density can be written as Seifert1 ; Lipowsky ; Seifert2

[TABLE]

where $\lambda_{i}$ is the Lagrange coefficient. The Lagrange equations \frac{\partial\Pi_{i}}{\partial\phi_{i}}-\frac{d}{ds}\big{(}\frac{\partial\Pi_{i}}{\partial\dot{\phi}_{i}}\big{)}=0 and \frac{\partial\Pi_{i}}{\partial x_{i}}-\frac{d}{ds}\big{(}\frac{\partial\Pi_{i}}{\partial\dot{x}_{i}}\big{)}=0 yield

[TABLE]

The Hamiltonian function is $H=-\Pi_{i}+\dot{\phi}_{i}\frac{\partial\Pi_{i}}{\partial\dot{\phi}_{i}}+\dot{x}_{i}\frac{\partial\Pi_{i}}{\partial\dot{x}_{i}}$. It has been proved that $H\equiv 0$ Julicher and it gives

[TABLE]

It also has proved that Eqs. 7, 10, 11 and 12 are identical to the general shape equation 3 Julicher . Note that $\gamma_{i}$ only appears in Eq. 12, so we can ignore Eq. 12 and only solve the left three equations because it can always be satisfied by choosing a suitable $\gamma_{i}$. This method will be used in following text.

In a multicellular system without any symmetry, whether the general shape equation 3 or the Lagrange equations 7, 10 and 11 are difficult to be solved analytically. In our work we suppose the system is composed by equal cells which have the same physical parameters: $\kappa\equiv\kappa_{i}$, $\omega\equiv\omega_{ij}$, $V\equiv V_{i}$ and the cell length $L\equiv\oint ds_{i}$. We will try to find analytical solutions with high symmetry for these Lagrange equations and use the finite element simulation method to study the complicate system and do stability analysis.

II.2 Analytical method without volume constraint

When cells adhere together and form periodic structures, there are many kinds of possible forms. Three typical forms are shown in Fig. 2. For simplicity, in each form we think the cells are equal and there is no volume constraint for each cell. In this case, we can make an ansatz that all of the adhesion parts are straight lines. (This assumption will be confirmed by finite element calculation in latter text). In these periodic structures, we can find the basic unit as shown in Fig. 3, which can be used to describe their equilibrium shapes.

In Fig. 3, due to the symmetry, $S_{1}$ and $S_{2}$ parts are equal, thus we can only choose the $S_{1}$ and $S_{3}$ to study. The Lagrange density in the region of $0\leq s\leq s_{0}$ for $S_{1}$ part is

[TABLE]

In the adhesion region $s_{0}\leq s\leq s_{0}+B/2$ , it is

[TABLE]

where $\kappa_{3}=\kappa_{1}+\kappa_{2}$. The shape equations for the $S_{1}$ part in $0\leq s\leq s_{0}$ can be written as

[TABLE]

For the $S_{3}$ part in $s_{0}\leq s\leq s_{0}+B/2$, there are

[TABLE]

At the initial contact point $s=s_{0}$, the boundary conditions in Eqs. 4, 5 and 6 are changed to

[TABLE]

Considering that all cells are equal, we have $\kappa_{1}=\kappa_{2}=\kappa_{3}/2=\kappa$. Note that the $S_{1}$ and $S_{2}$ parts should be equal in a periodic system, and that the rotation directions of them are opposite in the adhesion region, at the point $s=s_{0}$, there are

[TABLE]

The above conditions make sure that Eqs. 20 and 21 can be satisfied, and Eq. 19 is reduced to

[TABLE]

Due to each adhesion part is a straight line, $\phi_{3}$ should be independent to $s$ and we choose $\phi_{3}\equiv\phi_{0}=\pi/n$. One can find that, $n=2,~{}4$ and $6$ are corresponding to the shapes in Fig. 2(d), (e) and (f), respectively. In order to obtain dimensionless results, we choose $\kappa=1$ and fix the girth of each cell $\oint ds=2\pi$ Zhou1 .

First, we study the structures without volume constraint for each cell, which means the pressure difference $\Delta p_{i}=0$. (The non zero pressure case will be discussed in latter text). For convenient, in the following text we choose $\lambda_{1}=\lambda$, $\phi_{1}=\phi$, and $x_{1}=x$. The shape equations for the $S_{1}$ part can be reduced to

[TABLE]

The first integral is

[TABLE]

where $C$ is an integral constant. Consequently, we have

[TABLE]

Choosing the initial conditions $x(0)=0$, $y(0)=0$ and $\phi(0)=0$ and using the Eq. 7, the coordinates at the point $s=s_{0}$ are

[TABLE]

where $\text{FN}[x,y]$ and $\text{SN}[x,y]$ are the first and second incomplete elliptical integral, respectively. The length constraint $L=\oint ds=2\pi$ gives the following equation

[TABLE]

It yields

[TABLE]

Making use of Eq. 25, condition 23 is changed as

[TABLE]

Substituting the above equation into Eq. 30, we get

[TABLE]

The diameters of the inscribed circle and circumscribed circle of each cell in Fig. 2 are

[TABLE]

The dimensionless total energy for each cell is

[TABLE]

Substituting Eqs. 31 and 32 into the above equation, we obtain

[TABLE]

In the above equation, if $\omega$ and $\phi_{0}$ are known, there is only an unknown constant $C$. Defining $\chi=1/C$, the equilibrium shapes satisfy

[TABLE]

We find this equation yields $\chi\rightarrow 0$ ($C\rightarrow\pm\infty$). Then substituting Eq. 31 into Eq. 25, we have

[TABLE]

This result indicates that the optimal shapes of the free parts in Fig. 2 are circular arcs with the same radius $r_{0}=\sqrt{1/\omega}$ and the total energy for each cell $E=2\pi\sqrt{\omega}-\pi\omega$. Then the equilibrium distance between two cell is

[TABLE]

Note that the total length of each cell is fixed to $2\pi$, it needs the equilibrium shapes satisfy $r_{0}<1$. It yields $\omega>1$, that is the basic condition to form the cell-cell adhesion structures. If the adhesion occurs between a cell and a rigid plane, it needs $\omega>0.5$ Seifert1 ; Zhou2 .

II.3 Simulation method and the software

To study the complicate adhesion system, we use the software Surface Evolver Brakke which has been generally applied to simulate the equilibrium shapes of single vesicle Yan ; Zhou3 , cell adhesion system Ziherl , random foam Kraynik2 and droplet adhesion structure Lv . First, it needs to build an initial geometric model for the structure. Here we give an example model for three cells’ adhered together in the Supplementary Material. Next, we need to define the energy functional of the system. For one cell it is

[TABLE]

where $\alpha$ is a constraint and we set $\alpha=1/2$. Also, we fix $L=\oint ds_{i}=2\pi$. Defining the reduced volume $v_{i}=V_{i}/\pi=\int\int d\sigma_{i}/\pi$ we have $V_{i}=\int\int d\sigma_{i}=\pi v_{i}$ and $0<v_{i}\leq 1$ due to $L=2\pi$. Then we can fix the volume $V_{i}$ by setting the value of $v_{i}$. In this software the adhering potential $\omega_{ij}$ is named “tension” of the edges. After setting these parameters: $\alpha,L,\upsilon_{i}$ and $\omega_{i}$, we can use the Surface Evolver to find the equilibrium shape. In this process the perturbation method is used to find the optimal deformation direction of the system step by step and finally induces the shape converges to the equilibrium state with lower energy. More details of this can be obtained from Brakk’s home page Brakke . When it reaches the equilibrium state, the corresponding $B_{ij}$, $\gamma_{i}$, $\Delta p_{i}$ and the energy $E_{i}$ can also be obtained. The above method is used to find the solution of Eq. 3 numerically.

Although $\delta\Omega_{t}=0$ yields the equilibrium shapes equation Eq. 3, we cannot make sure each solution of Eq. 3 is stable because stable shapes need $\delta^{2}\Omega_{t}>0$. In the Surface Evolver, it can calculate the so called Hessian, which is the matrix of the second order differential coefficient of the total energy. At an equilibrium point, if the Hessian is positively defined, this means the point is a strict local minimum. Here we can give a simple explain about this method. In the finite elemental method, a continuous curve is divided into many straight lines. Then the elastic energy functional can be written as $\Omega(x_{i})$, where $x_{i}$ is the coordinates of the vertices of two neighbor straight lines. At $x_{i}=x_{i0}$, we expand the energy functional $\Omega(x_{i0}+\Delta x)=\Omega(x_{i0})+a_{i}\Delta x_{i}+b_{i}(\Delta x_{i})^{2}+\cdots$, where $a_{i}=\frac{\partial\Omega}{\partial x_{i}}\mid_{x_{i}=x_{i0}}$ and $b_{i}=\frac{\partial^{2}\Omega}{\partial x_{i}^{2}}\mid_{x_{i}=x_{i0}}$. If the shape reaches an equilibrium state, there is $a_{i}=0$. Further, a stable shape also needs $b_{i}>0$. Acutely, $x_{i}$ should be a vector and $b_{i}$ should be a matrix which is named the Hessian in the Surface Evolver. More details of the Hessian matrix can be obtained from the introduction file of this software Brakke . After a longtime simulation, the shape will converge to an equilibrium state. Then using the Hessian order we can obtain wether the Hessian is positively defined or not. If it is, we can say the shape is stable. But if it is not, the shape possibly is not stable and needs a longtime simulation again.

III Analytical Results

III.1 The periodic bead-like shape without volume constraint

When $\phi_{0}=\pi/2$, we get the bead-like shape in Fig. 2(a). For the optimal shape we define the elastic coefficient

[TABLE]

The corresponding Young’s modulus is

[TABLE]

where $H_{0}=H(\chi\rightarrow 0)$. Fig. 4(a) shows the relationship between Young’s modulus and $\omega$, which indicates that $Y$ is nearly linear with $\omega$. Furthermore, we can define the Poisson’s ratio

[TABLE]

It is a constant and we hope this results can be tested by the future experiments.

Now we consider that the equilibrium bead-like structure is under the action of a couple of forces along the length directions and the structure is in balance. Fig. 2(d) shows the diagram of one cell constrained by a couple of forces. Then the force is

[TABLE]

and the strain is

[TABLE]

Where $F=F(\chi)$ and $\xi=\xi(\chi)$ only depend on $\chi$. So we obtain the relationship between $F$ and $\xi$, which is shown in Fig. 4(b) for different $\omega$. However, the adhesion structure will be failure when the adhesion length $B$ reduces to zero. Letting $B=0$ in Eq. 32, we get $\chi=\chi_{m}$. Substituting it into Equation 45, we obtain the maximal strain $\xi_{m}=\xi(\chi_{m})$. Fig. 4(c) shows the relationship between $\xi_{m}$ and $\omega$. Also, we can get the pull-off force $F=F(\chi_{m})$ for different $\omega$ as shown in Fig. 4(d).

III.2 The square shape and the hexagonal shape without volume constraint

When $\phi_{0}=\pi/4$, we obtain the square shape in Fig. 2(e). For the optimal structure, we define the elastic coefficient

[TABLE]

If there are small perturbations $\Delta R$ along the square edges, the corresponding force is $\Delta F=K_{S}\Delta R$. The pressure on each edge of the square is $\Delta P=\Delta F/R_{0}=K_{S}\Delta R/R_{0}$. The valid area (the area for the circumscribed square) for the cell is $A_{0}=R_{0}^{2}$ and the area perturbation $\Delta A=2R_{0}\Delta R$. So, the bulk modulus (actually the area modulus in 2D case) is

[TABLE]

When $\phi_{0}=\pi/6$, using the similar method, we obtain the elastic coefficient as

[TABLE]

The bulk modulus for the optimal shape in Fig. 2(f) is

[TABLE]

We show the bulk modulus in Fig. 5(a). Apparently, with the similar $\omega$, the hexagonal structure has a bigger bulk modulus than the square structure.

As shown in Fig. 2(e), two couples of forces are supposed to act on this square structure and keep it in balance. The pressure on each edge of the circumscribed square derived from the out side forces is

[TABLE]

Correspondingly, for the hexagonal structure in Fig. 2(f), the pressure on each edges of the circumscribed hexagon is

[TABLE]

The volume strain for the two kinds of shapes induced by the force is $\sigma=R^{2}/R_{0}^{2}-1$. Fig. 5(b) shows the relationship between pressure $P$ and the strain $\sigma$ for the square structure (dashed lines) and the hexagonal structure (solid lines). Let $B=0$, we obtain the maximal strain shown in Fig. 5(c), which yields the adhesion failure boundaries (the two black lines) in Fig. 5(b). Moreover, we show the maximal strain and failure pressure for each structure under different adhesion potential in Fig. 5(d), which can be taken as the failure strength for these structures.

IV numerical results

IV.1 The multicellular bead-like structures

First, the multicellular bead-like structure are studied. In Fig. 6, each structure is composed by 10 cells with $\omega=8$. Gradually changing $v$, the symmetry breaking is found. When $v<0.75$, the formation is a straight line and has $D_{2}$ symmetry. When $v>0.75$, the structure will bend and presents periodic formation. Following the increase of $v$, the amplitude of the structure will increase and the period length tends to decrease. It means the cell number contained in one period is negatively related to $v$. Similar symmetry breaking induced by changing adhesion potential ware found between two red blood cells Ziherl . For fixed $\omega$, we find there is an optimal $v$, at which the energy reaches its minimum. As shown in Fig. 7, the total energy has the minimal value nearby $v=0.6$. A theoretical explain will be shown in latter text.

Fig. 8 shows the deformation induced by the change of adhesion potential at $v=0.94$. We can see that, with the increase of $\omega$, the amplitude and the cell number contained in one period will decrease. These waved 2D structures in Fig. 6 and Fig. 8 imply that the bead-like 3D shapes possibly present curved formations as well as helical structures. Some kinds of bacteria and algaes present bead-like structures, such as streptococcus.

IV.2 The dense stacking systems

Fig. 9 depicts two kinds of formations which are assumed to investigate the square adhesion systems in which each cell adheres to its four neighbors and the system can be extended to the whole 2D space. For Fig. 9(a), the shape appears in the region $0.94<v<1$ and always keeps the $D_{4}$ symmetry. When $v<0.94$, the free parts of the central cell will adhere to the other four neighbor cells which do not adhere to it before. As to Fig. 9(b), the stable region is $0.96<v<1$ and it also has the $D_{4}$ symmetry in this region. When $v<0.96$, the four cells in the corner of the out square will adhere to their nearest neighbors in the middle layer.

Fig. 10 shows two kinds of structures composed by 7 cells. The shape in Fig. 10(a) nearly has $D_{6}$ symmetry. Gradually decreasing $v$, the structure changes to the shape in Fig. 10(b) at $v=0.96$ with $D_{1}$ symmetry. When $v<0.92$, the two outer cells in Fig. 10(b) will lose the adhesion to the central cell due to the $m_{1}$ and $m_{2}$ adhesion parts reduce to zero. Fig. 11 shows two kinds of structures composed by 19 cells. The shape in Fig. 11(a) nearly has $D_{6}$ symmetry. But following the decrease of $v$, the non-symmetrical configuration will become manifest and the shape will change to Fig. 11(b). When $v<0.98$, our simulation indicates that there are many different shapes which have very similar total energy and the shape will change from one to one following the increase of simulation time. We think they are degenerate energy shapes and the energy barriers between them are very low. The degeneracy will increase with the decrease of $v$. When $v<0.92$, the two outer cells in Fig. 11(b) will lose the adhesion to the middle layer cells due to the $t_{1}$ and $t_{2}$ adhesion parts reduce to zero.

Our simulation also reveals that, when old adhesion parts become separate, new adhesion parts will probably appear somewhere. Unfortunately, at present time it is very difficult to track the evolvement of adhesion because the adhering or departing process between two cells can not be achieved freely in simulation. One possible way is to construct every possible formation and simulate their evolution processes. So, we can obtain the whole evolutionary path by connecting each independent process. However, this method is extraordinarily time-consuming when a stacking system contains plentiful cells.

IV.3 Periodic Shapes

In the former section, periodic shapes without volume constraint is studied analytically. If we consider the volume constraint, it is very difficult to obtain analytical results. Numerical simulation is a feasible way to find stable shapes. But at first, we tested our former analytical results with $\Delta p_{i}=0$ by the Surface Evolver. Our simulation indicates that, without volume constraint, all of the adhesion parts in the periodic structures in Fig. 2 are straight lines and our former analytical results with $\Delta p_{i}=0$ are reliable.

Further, three kinds of periodic structures with volume constraint for each cell are investigated. In the Surface Evolver, obtaining a stable periodic shape needs a long time searching. Fig. 12(a) shows a stable bead-like structure, in which a period contains one cell. To attain this shape, we need to scan the centre-to-centre spacing $R$ between two neighbor periods (note that a period can contain two or more cells) and find out the optimal distance $R$ at which one period has the lowest total energy. An example is shown in Fig. 13. Particularly, lots of simulations indicate that the optimal distance is close to the following analytic value. For a given $v$, supposing each cell is composed by two half of circles and two straight lines like the shape in Fig. 2(d), the distance between two neighbor cells is

[TABLE]

and the total energy for each cell is

[TABLE]

For example, when $v=0.8$, we have $R_{s}=1.10$ and $E_{t}=0.06$, which are close to the simulation values in Fig. 13. Besides the staircase shape in Fig. 12(a), it also reveals that there is another stable bead-like shape as shown in Fig. 12(c), in which a period contains two similar cells. Between these two stable shapes, it is the critical structure shown in Fig. 12(b) that one cell is composed by two half of circles and two straight lines. The centre-to-centre spacing between two neighbor cells satisfies $R=R_{s}$. But this critical shape is unstable. Fig. 13 shows the energy-distance relation obtained by simulation, which indicates there is an energy barrier between the two stable structures. The critical shape is on the unstable saddle point and the energy change in phase transformation is discontinuous. Therefore, it is possibly a second-order phase transformation between the two stable shapes and can be induced by the change of distance. Simulation also indicates this phase transformation can be found in the whole volume region ($0<v<1$). Actually, similar results have been reported Ziherl to explain the erythrocyte stacking structures.

If $\omega$ is fixed, the energy in Equation 53 has the minimum at a suitable volume. Let $dE_{t}/dv=0$, we get

[TABLE]

Fig. 14 shows the curve of $\omega$ vs $v$. Our simulation also agrees with this theoretic result. For example, Fig. 3 has shown that the total energy for the 10 cells’ bead-like structure has the minimal value nearby $v=0.6$ with $\omega=8$. Choosing $\omega=8$, Equation 54 gives $v\approx 0.58$.

When comparing the periodic shapes in Fig. 12 and the multicellular adhesion structures in Fig. 6 and Fig. 8, there are two evident differences. First, there are $D_{2}$ symmetrical structures for multicellular adhesion system, such as Fig. 4(a), (b) and (c). But we didn’t find stable $D_{2}$ symmetrical formations for periodic bead-like shapes. Second, the staircase configuration in Fig. 12(a) only appeals in the periodic structures and there is no similar stable shapes for the multicellular shapes.

Very specially, if Equation 39 is satisfied, the critical shape which is composed by two half of circles and two straight lines will be an analytical solution. Then $R_{0}=R_{s}$ coincidentally yields the same equation to Equation 54. Therefor, Equation 54 gives an analytical global optimal solution for periodic bead-like shapes. When Equation 54 is satisfied, our simulation reveals that the energy barrier is vanished and the critical shape is a stable middle phase between the staircase structure and the zig-zag structure. Therefor, the energy barrier only occurs when Equation 54 can not be satisfied.

Besides the bead-like shapes, two other periodic formations are investigated. Fig. 15(a) shows a stable rhombus structure which has two periods on two upright directions. The Surface Evolver code for this model is shown in the Supplementary Material. Our simulation indicates that the stable shape needs the periodic length on two upright directions are equal. It also reveals that the optimal distance is close to the following analytic value. For a given $v$, supposing that each free part of the rhombus shape is a quarter of circle with the same radii and the adhesion parts are equal straight lines, the distance between two neighbor cells is

[TABLE]

It gives $R_{s}=1.773$ for $\omega=4$. Fig. 15(a) is very close to this critical state. If $R<R_{s}$, the shape Fig. 15(a) will change. But what kind of structure will form is still unknown. If we gradually decrease $v$, the central cell’s two untouched neighbors will adhere to it when $v<0.92$. Then new phase will form and in which each cell adheres to it’s six neighbors, such as Fig. 15(b). From Fig. 15(a) to (b) the $D_{2}$ symmetry remains unchanged. For Fig. 15(b), the adhesion regions clearly are no longer the straight lines.

When structures are periodic in three directions, the calculation will be difficult because it needs to search the optimal distance in three directions. Our simulation indicates optimal shape possibly has unequal periods on different directions. Moreover, complex formations will appear, in which one periodic unit will contain many cells. It leads the problem to be more intricate. Therefor, the study of complex stacking system is still in challenging.

V Conclusions

We have studied the cell stacking system due to adhesion in 2D case. Compared with the random foam model, our model is more fit the reality of biotissues. We derive the analytical results for three kinds of periodic formations without volume constraint. The corresponding mechanical parameters, such as the Young s modulus, bulk modulus and failure strength, are obtained. We provide a potential method to connect the microstructure with the macro-mechanical parameters of biotissues. We also find that the increase of cellular volume will induce the $D_{2}$ symmetrical bead-like shape changing to squiggly shapes. These waved shapes imply there may be stable helical formations in 3D case. Our simulation reveals that there are plentiful degenerate energy shapes when many cells adhere together. The adhesion systems are flexible to transfer between these degenerate states. For periodic bead-like shape and foursquare shape with fixed volume, the lowest energy formations are nearby the critical states. Especially, there is a globally optimal solution for periodic bead-like shapes. We think the above results are helpful to understand the physical mechanism for the formation of biotissues. But it needs a more complicate investigation in 3D case, which will be our future work.

Acknowledgements.

We would like to thank professor Ken Brakke for his kind help on the use of the Surface Evolver. This work is supported by the National Natural Science Foundation of China Grants 11304383, 11304241 and 11374237.