Force Estimation during Cell Migration Using Mathematical Modelling

Cell migration is essential for physiological, pathological and biomedical processes such as, in embryogenesis, wound healing, immune response, cancer metastasis, tumour invasion and inflammation. In light of this, quantifying mechanical properties during the process of cell migration is of great interest in experimental sciences, yet few theoretical approaches in this direction have been studied. In this work, we propose a theoretical and computational approach based on the optimal control of geometric partial differential equations to estimate cell membrane forces associated with cell polarisation during migration. Specifically, cell membrane forces are inferred or estimated by fitting a mathematical model to a sequence of images, allowing us to capture dynamics of the cell migration. Our approach offers a robust and accurate framework to compute geometric mechanical membrane forces associated with cell polarisation during migration and also yields geometric information of independent interest, we illustrate one such example that involves quantifying cell proliferation levels which are associated with cell division, cell fusion or cell death.


Introduction
Cell migration is a fundamental cellular process that is essential to life and is linked to many important physiological and pathological events such as the immune response, wound healing, tissue differentiation, embryogenesis, and tumour invasion [1][2][3][4][5][6][7].
During migration, mechanical processes play a pivotal role, for example cellular biomechanics direct its physical behaviour, as well as its cellular functions in the biological context of health and disease [8][9][10]. Cells also physically interact with their extracellular environments via mechanical forces, for example, cell division, apoptosis, bleb and mitosis [11][12][13]. The strength of the forces varies as the sensitivity of a cell evolves with surrounding biomechanical and biochemical stimulus [11].
A key determinant of cellular biomechanics is the actin cytoskeleton [8]. It contains dynamic actin architectures that continuously re-arrange and turnover. The cytoskeletal forces are exerted on a plasma membrane, which define and insure the stability of the interior of the cell [14]. At the leading edge, a protrusion force is generated by actin architectures. Membrane tension balances these locally imposed forces and ensures rear retraction [10]. The characteristic time scale is short, often sub-seconds. The measured forces suggest they may range from Pico-Newtons (pN) to Micro-Newtons (µN) [10,14]. The cellular force generation intertwines with many other processes, forming a complex system. In addition, a noticeable change in a single cell behaviour may lead to a significant event on a tissue scale, so unravelling the mutual interplay between physical interactions such as protrusion and retraction forces are essential to understanding cell dynamics [10].
According to the research by Lieber et al. in 2013 [10] and Barbieri et al. in 2021 [8], there is very little understanding and limited ways to quantify cellular forces. The former claims a hypothesis on how membrane tension is set and regulated by cells, but states there is very little evidence to either support or disprove it; the latter describes challenges in quantification due to technical constraints.
Traditionally, our understanding of cell dynamics often comes from visual inspection using high-throughput, high-resolution microscopy and related imaging techniques [8,15]. Phase-contrast microscopes utilise partially coherent illumination to extract quantitative phase data [16,17]. Interferometric-based techniques make use of the Fourier decomposition [18,19]. Other alternative techniques include optical coherence tomography [20] and digital holographic microscopy [15,21]. Additionally to imaging, Simson et al. in [9] reports an interferometric technique to measure bending modulus, membrane tension and adhesion energy. A mechanic-optical biosensor is described in [22] to sense local cell adhesive forces. In [10] the authors discuss about membrane fluctuations and [23] summarises soft polymers that are typically used to measure cellular forces.
In this study, we propose an alternative theoretical and computational approach whereby, instead of measuring physical quantities in experiments, we describe the underlying rules (and often hypotheses) using mathematical equations, thereby obtaining a model for a migrating cell which incorporates certain assumptions on the physics underlying migration. There have been a number of studies carried out using simulations of such models to model cell migration, e.g., in [24] and subsequent related works, a phase-field model for keratocyte migrations is developed, in [25], some quantitative predictions are derived on how adhesion geometry and stiffness change cell behaviour. In this paper, we approach the problem of membrane force estimation during cell migration as the problem of computing forces such that we fit an established model of cell migration [26,27] to microscopy data that provides the cell membrane position at a series of time points. Specifically, we use the frames of imaging data to extract the position of the cell membrane at a series of times and use this data in an optimal control model as our target positions for the position of the cells under our mathematical model at the corresponding times. The control which is computed to minimise the difference between the cell positions generated computationally and the data corresponds to the protrusive force active at the cell membrane. This approach, i.e., the optimal control of phase field models albeit in a different context has received recent interest e.g., [28]. Computational simulations also help to build devices which can then be used to directly measure cellular properties, such as the microsystems summarised in [11].
The major novelties of the present work with respect to [26,27] are threefold: firstly, an application of the optimal control approach to time series data from real experiments, secondly, parameterisation of the model allowing the control problem to be interpreted as an approach for the estimation of forces during cell migration and finally, the application of the approach to the biologically important problem of quantifying cell division (or apoptosis) rates in a population of cells. These novelties, transform the approach from something that is of mainly theoretical interest to something that is of considerably utility to biological practitioners. This paper is organised as follows. In Section 2, we describe our theoretical and computational modelling approach. We take experimental observations of cell migration from three different cell types: keratocyte in [29], epithelial bladder cancer cell from the T24 cell line in [30], and epithelial kidney cancer cell from the MDCK cell line in [31]. Using our theoretical model, we re-create the corresponding computational cells and compute the predicted membrane forces. In Section 3, we present our results. We conclude our results in Section 5.

Mathematical Model for Membrane Force Estimation
We consider a volume-constrained Allen-Cahn equation with forcing as the "forward model" in an optimal control approach to whole cell tracking. The model arises from considering a simplified force balance on the cell membrane, further details on the model, the corresponding sharp interface formulation and its physical justification are provided in [26]. The volume constrained Allen-Cahn equation with forcing, a diffuse interface approximation to forced mean curvature flow, is stated as follows In the above model, φ( x, t) is a phase field variable whose zero level set corresponds to the cell membrane; , satisfying 0 < 1 is a parameter governing the interfacial width of the diffuse interface, G(φ) = 1 4 (1 − φ 2 ) 2 is a double well potential which has minima at ±1 and λ is a time-dependent constraint on the mass of the cell that models a volume constraint [32]. In practice our constraint differs from conservation of mass since the target data may have differing 'mass' (effectively cell area) from image to image. The initial condition φ 0 ( x) is taken as the initial image from the experimental observations. ν Ω is the outward normal to ∂Ω.
Without loss of generality we assume the domain Ω = [0, L] 2 . Biologically, φ( x, t) can be viewed as a volume fraction (φ ≈ 1 in the cell bulk, φ ≈ −1 in the extracellular matrix, φ ≈ 0 on the cell membrane), as the thickness of cell membrane, τ is an effective friction due to the interaction with the extracellular matrix, δ is the surface tension and since we focus on the two-dimensional, rather than three-dimensional cases, we average assuming a constant cell height of 0.1 µm.
In our modelling framework, η( x, t) is a membrane force generated by the cell during migration. A positive η indicates a protrusive force that drives the cell forward, while a negative force corresponds to a retractive force that allows the cell to contract enabling the cell to dislocate from the substrate to move its body forward. We make the assumption that the cells move as a result of forces that are only exerted in a region close to the membrane which is biologically reasonable since forces leading to migration are primarily exerted in the the actin cortex which is a thin region close to the membrane [26,27].
The physical interpretation of the model variables in Equation (1) is given in Table 1 below.
In the above equation, the units are µm : micrometer, s : second, and pN : pico-newton. In Table 2 we present the values we use for the parameters in the model, together with the references that gave rise to these choices. From here onwards, we wish to work with a unit-free model, and to this end, we introduce the dimensionless variables where the characteristic scale F is related to the kinetic scaling.

Formulation and Approximation of the Optimal Control Problem and Biological Interpretation
The approach we consider can be thought of as finding a forcing termη(¯ x,t) in Equation (4) such that the model best-fits the images. The methodology was first proposed in [26]; it sought to findη(¯ x,t) that minimises the following objective functional where N obs is the number of images we wish to fit to, φ obs,i (¯ x) is a (phase-field) representation of the data we wish to fit to at time t i and θ > 0 is a regularisation parameter that is necessary for the optimal control problem to be well posed. We omit the details on the formulation and solution of the optimal control problem in this primarily applications focussed work, referring the interested reader to references [26,27]. The primary computational work in solving the optimal control problem lies in approximating the solution to Equation (1) which (typically) must be carried out a number of times. We refer the reader to [34][35][36] and references within for more details on phasefield models and their solution methods. We previously developed two approaches for the approximation of the optimal control problem, a finite element approach in [26] and an adaptive, parallelised finite difference approach in [27], which is based on geometric multi-grid methods. The computational cost in fitting to the multiple datasets that we consider below prompts us to use the more efficient approach of [27] in this work and we refer to [27] for further details on the space-time discretisation. For completeness, we briefly describe our approach for obtaining the results presented in this paper in the Appendix A.
We note that the computedη(¯ x,t) can be interpreted physically as the membrane force (protrusive or retractive) required such that the motion of the boundary of the simulated cells under the model best resembles that of the imaging data. The interpretation of this η(¯ x,t) we employ in the present work is that it corresponds to the estimate of the forces exerted on the membrane, e.g., F-Actin based protrusion or myosin based contraction which govern the motion of the cell.

Model Parameters for the Different Biological Datasets
In Table 3 we give the spatial and temporal settings associated with the three cell types to which we apply our model: keratocyte in [29], epithelial bladder cancer cell from the T24 cell line in [30], and epithelial kidney cancer cell from the MDCK cell line in [31].

Results
For each of the biological datasets considered, we treat the initial frame of the video as data used to generate the initial conditions for the model. The remaining frames in the dataset are the target data we seek to fit the model to. Further details on our approach to extracting a phase field representation of the cell from imaging data are given in [26,27]. Our approach gives us the computed cell positions together with the estimated force such that the motion of the cells under our model recreates the observed motion from the imaging data. We use Figure 1 to illustrate our model when applied to experimental data from T24 cell line [30]. There is an 8-min gap between two adjacent frames in this experiment, and we take frames 3 and 4, for example. Our discretisation yields 10 time steps between these two frames. The first row in Figure 1 has two adjacent frames from the experimental data, and the second row shows the initial computed cell outline (obtained from the previous computation covering frames 2 to 3) and solutions at the 10 time steps with the computed optimal force. The solution from the 10th time step would then be used as the initial shape to compute the next stretch between frames 4 and 5. The dark shadow in the background shows the target shape as the objective, which is the shape of the cell from frame 4. This process continues successively throughout the full dataset. Figure 1. The first row illustrates two adjacent frames from the T24 experimental data [30] that were taken 8 min apart. The second row shows the initial shape and computed solutions at 10 intermediate time steps accordingly. The dark shadow in the background shows the targeted shape as the objective, which is the shape of the cell from frame 4. Bars indicate 20 µm.
As an example in Figure 2, we present the first frame of our results on T24. We show the original image from T24 experiment [30] on the top-left; our segmentation of the shape of the T24 cancer cell on the top-right (this segmentation technique is a combination of Otsu and edge detection, we refer the reader to [27,37] for more details); on the bottom-left, we demonstrate the interfacial region of the cell, its centroid position, and we continually overlay the cell shapes as the cell migrates. On the bottom right, we show the exerted forces where we use colour coding (red as protrusion and blue as retraction) to illustrate the location and amount of forces exerted on the cellular interfacial region, representing the cell membrane In Table 4, we summarise both protrusion and retraction forces re-created during the simulations of keratocyte migration (shown in Figure 2 and video in Appendix A). Each experimental image serves as a starting position or a goal. Our estimated forces are evaluated between adjacent frames. For the keratocyte [29], the actual real-world time between frames is 20 s. In Table 4, we show the average cell membrane length from our simulation in µm, the accumulated forces, the number of reconstructed time steps (denoted by RTS throughout), and the percentage of cell membrane where protrusion or retraction forces are exerted. Table 4. Details of estimating the membrane forces and its evolution through cell morphology reconstruction of the keratocyte. C. M. is an abbreviation for cell membrane. We note here that the Lagrange multiplier λ(t) in (1) is in effect a global spatially constant volume constraint force which is not included in the totals stated above. Our results on the epithelial bladder cancer cell T24 [30] are shown in Table 5. The layout of the table and its corresponding video in Appendix A are very similar to the keratocyte simulation. Table 5. Details of estimating the membrane forces and its evolution through cell morphology reconstruction of the T24. C. M. is an abbreviation for cell membrane. We note here that the Lagrange multiplier λ(t) in (1) is in effect a global spatially constant volume constraint force which is not included in the totals stated above. The results of the epithelial kidney cancer cell MDCK [31] are shown in Table 6 and are presented in the similar manner, apart from an additional diagram shown in the corresponding video in Appendix A on the right-hand side. Table 6. Details of estimating the membrane forces and its evolution through cell morphology reconstruction of the MDCK. C. M. is an abbreviation for cell membrane. We note here that the Lagrange multiplier λ(t) in (1) is in effect a global spatially constant volume constraint force which is not included in the totals stated above. We note that making direct comparisons of our results with experimental studies is challenging since other forces relevant to migration, such as traction forces generated through interactions with the substrate, are neglected in our model and often membrane forces are not directly measured. However, our results are consistent with available experimental results which attempt to measure retractive and/or protrusive forces generated by migrating cells; these results conclude that the total force exerted by the cells is of the order of 10s of Nano-Newtons [38].

Geometric Quantities That Are of Biological Interest
Based on the computed phase field information, we can compute important and biologically relevant geometric information that is of significant interest to experimentalists. Through simple post-processing of the model outputs, we may obtain geometric quantities such as circularity, curvature and elastic energy of the membrane. We demonstrate one such post-processed quantity of interest by considering a dataset that consists of multiple cells undergoing division and demonstrating how our approach allows us to quantify cell proliferation rates, and also demonstrate geometrically, the cell division process. To proceed, we first state how to compute the Euler number of the cell membranes which in effect corresponds to the total number of cells present in the simulation. The Euler number, in two dimensions in the phase-field formulation is given by [39] Here φ is the computed phase-field function. As this number corresponds to the number of cells present, it is extremely useful as it gives us the means to automatically track events such as cell division and cell fusion during the process of cell migration. This approach is extremely valuable as it could automate an otherwise laborious task of counting division, fusion or death events and removes the need for genetic manipulations which would be required to highlight such events. We illustrate this in the video in Appendix A related to the kidney cancer cell MDCK [31] where a number of cell divisions occur during the video and these are tracked accurately by the Euler number of the computed phase field.
For completeness, we include the final frame of the result video on the kidney cancer cell MDCK in Figure 3. In this figure, the original image from the experimental observation is shown as the first image on the first row, with a red box highlighting our choice of three cells used in our simulation. The second image on the first row illustrates the segmentation of the corresponding cells in the first image on the first row. The first image on the second row shows the interfacial region of the selected cells and their centroid points. In this sub-figure, we continually overlay the cell shapes and portions as they migrate. We use red for protrusion and blue for retraction forces to identify the regions where they are exerted around the cellular interfacial region. The dark shadows in the background illustrate the targeted shapes our model is replicating, and the bar on the right-hand side shows the maximum and minimum amount of forcing that the colour coding is illustrating. The only image on the third row shows our Euler number from Equation (5) computed at each RTS. Within this data, the initial three cells are divided into six cells, and we use red circles to indicate the cell division events in this graph.

Discussion
Single or population cell migration is essential for many biological processes such as immune response, embryogenesis, gastrulation, wound repair, cancer metastasis, tumour invasion, inflammation and tissue homeostasis. However, aberrant or defects in cell migration lead to various abnormalities and life-threatening medical conditions [40][41][42]. Increasing our knowledge on cell migration can help abate the spread of highly malignant cancer cells, reduce the invasion of white cells in the inflammatory process, enhance the healing of wounds and reduce congenital defects in brain development that lead to mental disorders. While single-cell sequencing has accelerated breakthroughs in cancer research and transformed our understanding of tumour biology, leading to significant impacts for cancer treatments, understanding and quantifying membrane force generation associated with cell migration remains an open problem, in this study we have exploited our modelling approach of using optimal control of surface geometric partial differential equations posed on the cell membrane to predict and estimate biological quantities of interest, such as protrusion and retraction forces. Our approach is substantially different from currentstate-of-the-art modelling of such forces, it sets premises to study cell migration through complex multi-dimensional environments where force generation between the cell and the extracellular matrix is critical [40][41][42].

Conclusions
A number of recent studies such as [10] remark that small changes in mechanical forces generated from individual cells can lead to fundamental changes at tissue levels. However, as [8] indicates, it is technologically challenging to simply measure those forces during experiments, such as during the process of cell migration.
In this work, we estimate the forces exerted by migrating cells by computing 'optimal' forces such that a mathematical model for cell migration best fits observed imaging data. In this paper, we took experimental data of three different cell types: keratocyte in [29], epithelial bladder cancer cell from the T24 cell line in [30], and epithelial kidney cancer cell from the MDCK cell line in [31]. For each case, we demonstrate how we re-create the observed cell migration and summarise the protrusion and retraction forces generated under our model. We also note our approach is applicable to multiple cells and can be applied in three dimensions, given appropriate datasets in 3D. Our approach could also allow us to access biologically relevant quantities such as membrane length, circularity and curvatures. We demonstrate one example using the MDCK cell line dataset [31] in which a number of cell divisions occur during the evolution. Our approach deals robustly with this setting allowing accurate quantification of cell proliferation rates which is generally cumbersome if carried out manually. Moreover, we provide a means of tracking (automatically) the number of cells present which could be of practical use if one wishes to measure the rate of cell divisions, cell death or cell fusion.
Our proposed approach is amenable to further improvements and these include computing more accurate measurements of parameters such as friction force and surface tension as well as more refined modelling of migration itself. We note that it would be relatively straightforward to extend our simulations to three space dimensions, which would enable the recreation of more accurate cells and their environments [27] but a major challenge in this case is obtaining sufficiently high-resolution imaging data. As [43] states, cell-matrix adhesions and cytoskeletal organisation could be different in 2D and 3D measurements, and may alter key cell responses, including morphology, migration and proliferation. We also note that in principle this approach can be adapted to more complex models of cell migration and this is merely a proof-of-concept study illustrating the utility of our approach. we use colour coding to identify red as protrusion, and blue as retraction forces and the locations they are exerted on the cellular interfacial region, the dark shadows in the background illustrate the targeted shapes that our model is trying to replicate, the bar on the right-hand side shows the maximum and minimum amount of forcing that the colour coding is illustrating.; the only image on the third row: we show our Euler number from Equation (5) computed at each RTS and red circles indicate the events of cell division. Bars indicate 20 µm. • Our approach to present simulation results The mathematical model in Equation (1) and its dimensionless form (i.e., Equation (4)) exist continuously in both space and time in the closed domain Ω. We apply discretisation schemes to obtain a discrete system at each reconstructed time step (RTS). For example, the situation presented in Figure 1 has been discretised into 10 time steps. Solving such a system yields its simulation result at its corresponding reconstructed time step. We use φ i and η i to represent the discrete solutions on a grid point i where there are N grid points, and Ω L is the number of grid points on one axis of the square domain Ω. φ i describes the evolution of the cell shape and takes values either −1 or 1 in-or outside of the cell. When φ i is between −1 and 1, it illustrates the phase-field interfacial region and the cell shape. The phase-field variable undergoes a smooth transition between the values of 1 and −1.
The membrane corresponds to the zero level set of the phase field variable. In order to evaluate quantities on the membrane in our simulations results we need to specify the interfacial region from our simulation by considering a small interval around zero to be the interfacial region since our method does not explicitly track the interface itself. We define grid points i that has values of φ i between −0.02 and 0.02 as the cell membrane. We evaluate the average length of the cell membrane on the interval I, where I = t T I using the formula, ∑ I * RTS k=(I−1) * RTS+1 ∑ N i=1 L Ω L | − 0.02 < φ k i < 0.02 /RTS. η i is the forcing, and we are interested in these exerted closely around the cell membrane. The total protrusion force on the interval I is evaluated using the formula, ∑ I * RTS k=(I−1) * RTS+1 ∑ N i=1 η k i | − 0.02 < φ k i < 0.02 and η k i > 0 . The total retraction force on interval I is evaluated using the formula, ∑ I * RTS k=(I−1) * RTS+1 ∑ N i=1 η k i | − 0.02 < φ k i < 0.02 and η k i < 0 . The choice of RTS affects the stability of the solution methods, and we refer the reader to [34][35][36] for the computational details.