A sub-cellular viscoelastic model for cell population mechanics

Abstract

Understanding the biomechanical properties and the effect of biomechanical force on epithelial cells is key to understanding how epithelial cells form uniquely shaped structures in two or three-dimensional space. Nevertheless, with the limitations and challenges posed by biological experiments at this scale, it becomes advantageous to use mathematical and 'in silico' (computational) models as an alternate solution. This paper introduces a single-cell-based model representing the cross section of a typical tissue. Each cell in this model is an individual unit containing several sub-cellular elements, such as the elastic plasma membrane, enclosed viscoelastic elements that play the role of cytoskeleton, and the viscoelastic elements of the cell nucleus. The cell membrane is divided into segments where each segment (or point) incorporates the cell's interaction and communication with other cells and its environment. The model is capable of simulating how cells cooperate and contribute to the overall structure and function of a particular tissue; it mimics many aspects of cellular behavior such as cell growth, division, apoptosis and polarization. The model allows for investigation of the biomechanical properties of cells, cell-cell interactions, effect of environment on cellular clusters, and how individual cells work together and contribute to the structure and function of a particular tissue. To evaluate the current approach in modeling different topologies of growing tissues in distinct biochemical conditions of the surrounding media, we model several key cellular phenomena, namely monolayer cell culture, effects of adhesion intensity, growth of epithelial cell through interaction with extra-cellular matrix (ECM), effects of a gap in the ECM, tensegrity and tissue morphogenesis and formation of hollow epithelial acini. The proposed computational model enables one to isolate the effects of biomechanical properties of individual cells and the communication between cells and their microenvironment while simultaneously allowing for the formation of clusters or sheets of cells that act together as one complex tissue.

Figure 1. A comparison of existing models for cell morphology based on model realism and computational cost.

A = Advanced, S = Simple, N = None; L = Low, M = Moderate, H = High.

Figure 2. Cell structural model.

a) The perimeter of the cell and nucleus (i.e. their corresponding membranes) are initially discretized into nodes (points). The superscript indicates that the point is on the cell membrane and superscript represents a point on the nuclear membrane. If neither nor are specified, the given point can be assumed to lie on either the cell membrane or nucleus. For example, the membrane's point of cell represented by . Each line that connects two points (red, green and blue lines) refers to a Voigt subunit. The total force that acts on each point is and is calculated by Eq(1) b) Voigt subunit. A linear Kelvin-Voigt solid element, represented by a purely viscous element (a damper) and purely elastic element (a spring) connected in parallel. The force that is exerted on from this subunit is (Eq.(2)). is the spring constant and represents viscosity.

Figure 3. Inner cell structure and forces.

The mechanical properties of the cytoskeleton are modeled using Voigt subunits; the spring constants of the model are linear approximations to the elasticity of the inner cell. All springs can be considered subject to a damping force due to the viscosity of the cytoplasm, where linear dash-pots are used to approximate the viscosity of the cytoskeleton. In our model, the cytoskeleton is divided into uniformly radial distributed parts, each of which is replaced by a Voigt subunit radiating from the nucleus (blue subunits). Each subunit connects two points of the cell and nuclear membrane, which are located at a radial direction from the center of the nucleus. The model also contains Voigt subunits in the nucleus (red subunits), each of which connect two nuclear membrane points and in which equal to , This allows the nucleus to show more resistance to changes in its shape and volume due to exterior pressure. is the cytoskeletal force acting on and is calculated by Eq. (3). is the force acting on from the cytoskeleton and nuclear cytoskeleton and is calculated by Eq. (4).

Figure 4. Structure of cell membrane and cytoplasm.

a,b) To represent the viscoelasticty of the membrane and cortical cytoskeleton, two consecutive membrane points are connected with a Voigt subunit (green subunits); hence, the model includes Voigt subunits on the cell membrane and subunits on the nuclear membrane. The forces acting on each cell from membrane subunits is calculated by Eq. (5); as the figures show, each point is subject to two adjacent subunits. c) An osmotic pressure will act on the membrane. This internal pressure is involved in cell morphology and affects the driving force of cell movement –. Knowing the persistence lengths of micotubles, and the fact that they appear curved in the cell, it follows therefore, that this filament pushes the membrane outward . Therefore, a pressure field acting upon each point of the cell membrane, representing cytoplasmic pressure with an outward and perpendicular direction to the cell membrane can be defined as by Eq. (6).

Figure 5. Mitosis and involved forces.

a) Several bio-mechanical aspects of cell proliferation are included in our model. First, the cell area (or volume in 3D) is doubled. The axis of cell division is selected (dash-line), in a dividing unpolarized cell, this axis usually is perpendicular to the cell elongation direction in such a way as to split the cell into two approximately equal parts. In a partially polarized cell, however, the axis of cell division is orthogonal to the part of the cell membrane that is in contact with the ECM. Two new daughter nuclei are then placed orthogonal to the axis of cell division. After the selection of division axis, the model finds the nearest membrane point to this axis, i.e. point In mitosis there are two major mechanical forces, first in the Anaphase stage of mitosis, shortening the spindle fibers caused by the kinetochores separation, and the chromatids (daughter chromosomes) are pulled apart and begin moving to the cell poles . Second, a contractile ring is formed by contractile forces acting on the opposite sites of the cell boundary in the cytokinesis process . This results in the formation of a contractile furrow and causes division of the cell into two daughter cells. The cell points can therefore be divided into two groups, A and B, where group A consists of membrane points from to and nucleus points from to and the remaining points belong to group B. To model the first mechanical force, the points of the nucleus in A and B are pulled apart, in the orthogonal direction to the division axis with force (Eq.(8) ). During the nucleus separation, the contractile force, , acts on boundary points of A and B groups to model the second mechanical force(Eq. (9)). b) The main phases of cell growth and division. (I) Cell growth. To implement cell growth in the proposed model, the number of membrane points, i.e. the number of viscoelastic compartments, is allowed to increase. When we add two points on each the cell and nuclear membranes, four subunits are added to the system, with the parameters of these new subunits calculated from the average of the first neighbor's homogeneous subunit parameters. With the additional ‘growth’ point, the circumferential length of the membrane increases in proportion to . Hence, the rest volume i. e. the volume of the cell when it grows freely without any inner or outer constraint, must increase proportional to . Therefore, the rest length of radial springs is increased in proportion to . When the area (or volume) of cell doubles, the number of defining membrane points increases to , where is the number of membrane points on the initial cell. (II,III) Mitotic process: two types of forces act on points to divide the cell. Due to these forces the cell elongates and prepares for division. (IV) Two new daughter nuclei are then placed orthogonal to the axis of cell division. After the nucleus separates, i.e. the distance between the center of the mass of nuclear points exceeds a certain value, , the cell will divide into two daughter cells, i.e. the subunits which join the boundary points will be eliminated and will bind to a new first neighbor point in the same group with a new subunit. V) After division takes place, each daughter cell will only have points, and as a result it is possible to simultaneously add points to each cell. To add membrane points, two consecutive points in the membrane are found that have the longest distance and a new point is inserted between them, and this process is repeated until the number of cell points becomes . f) Adhesion of the two daughter cells.

Figure 6. Motility.

Cell crawling is generated by the interplay between three different processes, namely, protrusion, adhesion, and contraction. These processes cooperate in a spatially heterogeneous structure to generate a complex topology for cell motion while correlation and coordination between them has a significant role on the motility of the cell . To model these three-stage events, the cell is first polarized by categorizing the points into two groups, the anterior and posterior, where their cytoskeleton subunit parameters will change periodically in a coordinated fashion. In addition to modeling the adhesion with a substrate, the drag coefficient is used for each of the points which will change periodically in coordination with variation of subunit parameters. The method for the variation of the parameters is shown in Eq.(10).

Figure 7. Apoptosis.

The structure and morphology of apoptotic cells undergo dramatic changes, including detachment from the neighboring cells, collapse of the cytoskeleton, shrinkage of the cell volume and alternations in the cell surface. Apoptosis progresses quickly and its products are quickly removed. To model these events, cell adherens junctions with their neighboring cells and/or substrate are initially disassembled (top right). Then, the spring constant length of all subunits are reduced arbitrarily and the inner pressure is removed (bottom left), so the cell will collapse and the cell area is gradually reduced until it reaches a prescribed minimal value (bottom right). At this time, the cell is considered to be dead and will be removed from the system.

Figure 8. ECM.

In 3D culture, our model allows for the investigation of a cross section of the system and the cells that are immersed in ECM. In 2D, the contact region of the ECM and the cell is a line that surrounds the cell; therefore, an enclosed curve (or ring) can be used for the ECM that surrounds the cells (blue curve). The ECM is modeled using a chain of subunits connected in series, where each point connects two subunits (blue curve with red points). These subunits can interact with cell points in a manner similar to the interaction of two points of different cells. This chain is flexible and the number of its corresponding points can be increased or decreased.

Figure 9. Cell-cell interaction.

Each cell can interact with another cell and substrate in two methods, adhesion and/or repulsive forces due to elasticity, as shown in Eq. (11). a) In this model, all points located on the cell membrane serve as potential sites of cell-cell connections. Each two points from different cells or cell-substrates are connected via a Voigt subunit, once they are closer than a determined value of . The adhesion subunit parameters between two cells are a function of these parameters, for more detail see Eq. (12). b) The repulsive force acts as a short range force. It is a passive force resulting from the elastic interaction with neighboring cells and acts on each point of the cell, when the distance to the other cell points or substrate is less than . The magnitude of the repulsive force is a function of the distance of two surfaces (Eq. (13)) and its direction is perpendicular to the membrane, pointing inward to the inner cell.

Figure 10. The monolayer culture of cells and the effects of adhesions on tissue formation and morphology.

a) In this simulation, D = 0.0001*D0 i.e. low intensity of adhesions. We began from two cells and allowed them to reproduce freely, subject to conditions that are suitable for proliferation without the occurrence of apoptosis. Results show a filled circular culture and fast proliferation. n represents the dimensionless elapsed time. b) D = D0 i.e. high intensity of adhesions. The process began from two cells which were allowed to reproduce freely. The results show dendrite morphology for the culture which can be seen frequently in tumors. It is conceivable that this special morphology is due to strong adhesion between cancer cells and the ECM. The process has a slower proliferation rate than part a. The epithelial cells in a monolayer appear as polygonal cells. It also can be seen that the average number of neighbors for any cell is 6, regardless of the value of the drag coefficient.

Figure 11. Effect of adhesion intensity on monolayer cell culture properties.

The value of for cultures with 120 cells, compared against the adhesion intensity. This graph shows a linear relation between and the adhesion intensity, which in turn suggests that the cellular shapes are almost circular at low drag coefficients, and diverge from being circular as the drag increases.

Figure 12. Growth of epithelial cell interacting with the ECM.

In this simulation a cross-sectional perspective of cell culture can be seen. Therefore, the ECM is a line of points to which the cells adhere. These points are dynamic and can change if needed. n represents the dimensionless elapsed time. a) Cells plated on a layer of surface culture. As the cell proliferates, a stable, uniform monolayer will be constructed. The axis of division is perpendicular to the ECM, likely related to polarized cells. If a cell detaches from the ECM due to the loss of polarity, it will activate the apoptosis pathway. b) A hole exists in the ECM, where a cell is located for attachment and proliferation. After polarization, the cell starts to proliferate and create a stable, lumen-containing cyst, lined by a single layer of epithelial cells. As it can be seen, the ECM is deformed a bit due to the dynamic interaction between the ECM and cells during the growth process. c) Shows an inverted cyst. A circular ECM is located in a suspended culture, to which a cell is attached and polarized. Upon completion of proliferation, cells surround the entire surface of the ECM and create inverted cysts, with matrix deposited on the inside of the cyst. If the process is allowed to continue, the cyst will grow further and become greater, which corresponds to a bigger ECM. This is because the volume of ECM in our model can freely increase.

Figure 13. Effects of a gap in ECM surface.

ECM is considered to be rigid and not affected by cells; however, cells still adhere to ECM and are polarized. a) Figure shows the final results for various gap sizes. If the gap width is denoted with and the radius of a free epithelial cell with , then it can be seen that cells cannot line the gap for . For the first cell which meets the gap will enter it, although due to the pressure of the walls it would not be able to continue its growth and division, so it fills the entry and blocks the gap. Other cells pass over the gap and again create a linear monolayer. For cells cannot ignore the gap and penetrate it. They continue their proliferation into the gap; however, when they reach the internal right corner, because of the limitation in space and the forerunning cells being subject to direction changing, the growth is stopped and the cells are entrapped in the gap. For the cells can enter the gap without any problem and line it. b) A few snapshots of the growth process when the gap is equal to . n represents the dimensionless elapsed time. During simulation, cells show differing behaviors at the corners. The growth rate of cells decreases at the internal corner and increases at the external corner. In addition, at external corners, due to the sudden decrease in contact area, cells detach from the ECM more easily in response to the pressure of neighboring cells. This, in turn, leads to loss of polarity and apoptosis.

Figure 14. Tensegrity and Tissue morphogenesis.

At this stage, an attempt is made to model the tensegrity hypothesis. To do so, a monolayer of epithelial cells is placed on the ECM. This causes the ECM to get thinner, by decreasing the drag coefficient at the center of the ECM (color gradient represents change in drag coefficient). The simulation is run with two densities of cells. The cells at this region have localized growth and motility, which drive the ECM downward. As a result, cells find space for more growth and proliferation. Therefore, they continue to drive the ECM further, which finally leads to the creation of a bud. This bud can be the first stage of tubulogenesis. a) Due to the high density of cells, they do not have enough space to grow and proliferate, so there is more organization and less deformation. b) Shows proliferation and less symmetry due to low cell density.

Figure 15. Formation of hollow epithelial acini.

This stage shows modeling snapshots of self-arrangement of individual eukaryotic cells into a stable hollow acinar structure. In this model a single cell (n = 1) undergoes several consecutive divisions and gives rise to a small cluster of cells containing two different populations (n = 100–217): the inner cells are entirely surrounded by other cells which do not have access to ECM, and the outer cells partially face the ECM. Further cell proliferation leads to the expansion of the whole cluster. During this stage (n = 150–217) some intercalation of outer cells to the interior for the preservation of the circular shape of the tumor can be seen. As the intensity of adhesion between cells increases, the process of intercalation becomes more difficult. After this stage, when the tumor reaches a certain age, the cell undergoes further differentiation of outer cells and results in their apical-basal orientation and self-arrangement into one layer of polarized epithelial cells of regular cubical shapes (n = 219) and the inner cells are then triggered by polarized cells to enter the process of cell apoptosis. Each cell that does not have access to the ECM and as a result does not get polarized will die (n = 219–225). This process leads to the creation of the inner lumen. Consequently, the proliferation of polarized cells is suppressed and the final structure stabilizes in the form of a hollow epithelial acinus(n = 250–450). Moreover, the processes of cell proliferation, polarization and apoptosis need to be well coordinated in order to maintain the hollow acinar structure in a stable manner; otherwise cell overgrowth may lead to intraductal carcinomas. This coordination shows that this process is very dependent on biochemical signaling between cells. The final shape of the tumor is very dependent on the viscosity of the ECM. If the viscosity of the ECM is high enough (i.e. drag coefficient is high in our model, due to this process trying to move during a minimal distance in energy space), the cell maintains its circular morphology. However, if the viscosity of the ECM is reduced, the cell deviates from the circular shape.