(C) PLOS One
This story was originally published by PLOS One and is unaltered.
. . . . . . . . . .



Derivation and simulation of a computational model of active cell populations: How overlap avoidance, deformability, cell-cell junctions and cytoskeletal forces affect alignment [1]

['Vivienne Leech', 'Department Of Mathematics', 'University College London', 'London', 'United Kingdom', 'Fiona N. Kenny', 'Randall Centre For Cell', 'Molecular Biophysics', 'King S College London', 'Stefania Marcotti']

Date: 2024-10

We build a mechanistic, agent-based model to describe the motion of individual cells interacting with neighbouring cells in 2D, where cells are approximated as ellipses. While typically elongated, real fibroblast shapes are of course more complex and include ruffles, and protruding and contracting lamellipodia [ 56 , 57 ]. We argue that, over the timescale we are interested in, ellipses are a rough, but appropriate approximation of real cell shapes. We also choose ellipses, because they allow for a straight-forward description of dynamic cell shape changes (see Sec 3). Further our modelling approach can be extended to more complex cell shapes, which will be the subject of future work. In this work we do not model cell divisions, however, they could be included in a straight-forward manner.

While alignment processes in active particles are relevant in many contexts, we will focus on the particular example of fibroblasts. Fibroblasts are cells in the connective tissue in animals and are responsible for making and remodelling the ECM. Fibroblast and ECM alignment is observed during various scarring pathologies. In [ 55 ] we have investigated the difference in alignment behaviour of fibroblasts in healthy tissue (normal dermal fibroblasts, NDFs) as compared to dermofibroblasts in certain scar tissue (keloid derived fibroblasts, KDFs). KDFs were found to show stronger alignment over larger length scales, Fig 1A and 1B . Further we found that KDFs show less tendency to crawl on top of each other and form aligned supracellular actin bundles via cell-cell junctions, spanning multiple cells. Using mathematical modelling, we found in [ 55 ] that the increase in overlap avoidance can explain the stronger alignment. In this work we will use the differences found in NDFs and KDFs to motivate further model extensions. However, the model ingredients are applicable to other cells and situations and hence the findings are relevant beyond fibroblasts.

2.2 Model derivation

We will show the derivation excluding cell deformability and cell-cell junctions. These will be considered in Sec 3 and Sec 4. More derivation details can be found in S1 Appendix, Sec 1 and a summary of model parameter names and meaning can be found in S1 Appendix Tab. 1. We consider N cells within the fixed domain , each with centroid position and orientation α i ∈ [0, 2π), i = 1, …, N. Each cell is described by an ellipse with semi-major axis a and semi-minor axis b as shown in Fig 1D. The cell’s area is given by A = abπ. In the absence of other cells, each cell self-propels with constant velocity w in direction e(α i ) = (cos(α i ), sin(α i ))T, where superscript T denotes the transpose. We assume that this self propulsion encompasses the mechanism of an individual cell moving in a directed way across the substrate. An alternative approach would be to model cell motion as a persistent random walk, as has been done in [58, 59].

We derive the governing equations using energy minimisation. An alternative model derivation based on force balance, which leads to the same governing equations, can be found in S1 Appendix, Sec 1. Focusing on one cell positioned at X(t) with orientation α(t) at time t, we parameterise the points inside the cell, x, by (1) where s encodes the distance from the centre of the cell X to the point x as a proportion of the distance from the centre of the cell to the boundary and θ encodes the angle parameter, see Fig 1C. The rotation matrix R(α) and the shape vector k(θ) are defined by We assume that at every time step Δt, the system minimizes a total energy E tot , which, for the base model, is the sum of contributions from friction E friction , from overlap avoidance E overlap and from self-propulsion E prop . All terms inside the integrals below represent the effect of each contribution on one point x inside the ellipse. We then obtain the total energy for one cell by integrating over whole (elliptic) cell area, i.e. with respect to s and θ. The chosen parametrisation given in (1) leads to the appearance of the area element abs (which accounts for the fact that closer to the cell centre one step in s-direction contributes less area) in the integrals in (2), (3) and (4). E friction models friction with the environment by comparing how much points have moved between time t and time t − Δt: (2) The overlap avoidance term E overlap is modelled by an energy potential V which includes overlap avoidance interactions with all other cells. The choice of V will be discussed below. (3) Finally, we want to model self-propulsion. One way to do this is to include it in the energy formulation by prescribing a force F acting on the cell. In the course of the derivation it is chosen to be F = wηe(α), i.e. acting in the direction of the orientation and proportional to the experienced friction η. This choice of F leads to a self-propulsion speed that is independent of friction (choosing F not proportional to η would only lead to a different definition of the non-dimensional quantities below). (4) The total energy is then given by summing (2), (3) and (4) (5) We obtain governing equations by minimising this energy in each time step. In other words, in each time step the cell can change its characteristics (position, orientation, shape) to decrease the energy. Calculation details can be found in S1 Appendix Sec 1. The main derivation steps for the base model are: 1. Differentiation with respect to X and α respectively (treating all other variables in the energy potential as constants). 2. Setting the derivative to zero and taking the limit Δt → 0. 3. Evaluation of the integrals. We then obtain the following differential equations for the motion of one cell (6a) (6b) The superscript ⊥ describes the left-turned normal vector. The two equations in (6) show how the position and orientations are influenced by the force and torque associated with V respectively. Note that we are working in a friction-dominated regime, which is why velocity and angular velocity (as opposed to acceleration and angular acceleration) are proportional to force and torque.

Choice of overlap potential V. The potential V describes the influence of overlap, where V > 0 describes overlap avoidance and V < 0 overlap preference. Many choices of V are possible: e.g. since cells might be thicker closer to the cell center, overlap closer to the cell center could be punished more than further away. However, due to the governing equations in (6) being formulated in terms of integrals of the gradient of the potential V, complicated shapes of V are computationally harder to evaluate, especially in the context of collective dynamics when this will need to be computed numerous times at each time step. We therefore choose V to be constant with value σ in regions of overlap and zero elsewhere: For two overlapping ellipses with domains and we define , where is the indicator function which equals 1 if and 0 otherwise. The strength of this potential is . If σ > 0, the cells experience repulsion in response to overlap, and if σ < 0, the cells experience attraction. In this work σ > 0. As a result of this choice of potential V, two cells only experience overlap avoidance upon overlapping with each other, hence we define as the set of indices of cells that overlap with the i-th cell.

Final base model. We non-dimensionalise the model using as reference time , as reference length and define r = a/b as the cell’s aspect ratio. The above choice of V allows to evaluate the integrals in (6) explicitly (for calculation details see S1 Appendix, Sec 1). The resulting equations can be formulated such that they depend only on the points of overlap between cells i and j, denoted by , where up to k = 4 points of overlap are possible. This is computationally advantageous since only the points of overlap need to be found, instead of areas of overlap which would be more computationally costly. In the following K ij = 1 or K ij = 2 denotes the number of overlap point pairs between cell i and cell j (having one or three points of overlap can be reduced to having zero or two points of overlap). The overlap points are ordered such that they traverse the boundary of the cell in an anti-clockwise direction. and (and and ) are chosen in such a way that the boundary segment of cell i between these pairs of points is contained in the domain of cell j, see Fig 1D. (7a) (7b) where the non-dimensional quantity ν is given by and can be interpreted as comparing the strength of repulsion in the presence of friction to the self-propulsion speed (see more interpretation in the results section below). These two governing equations are supplemented with initial conditions and boundary conditions. Throughout this work the domain is a square box with side length L and cells are initially placed randomly inside the box with a random orientation. Further, we use periodic boundary conditions.

[END]
---
[1] Url: https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1011879

Published and (C) by PLOS One
Content appears here under this condition or license: Creative Commons - Attribution BY 4.0.

via Magical.Fish Gopher News Feeds:
gopher://magical.fish/1/feeds/news/plosone/

You are viewing proxied material from magical.fish. The copyright of proxied material belongs to its original authors. Any comments or complaints in relation to proxied material should be directed to the original authors of the content concerned. Please see the disclaimer for more details.