Image based validation of dynamical models for cell reorientation

A key feature of directed cell movement is the ability of cells to reorient quickly in response to changes in the direction of an extracellular stimulus. Mathematical models have suggested quite different regulatory mechanisms to explain reorientation, raising the question of how we can validate these models in a rigorous way. In this study, we fit three reaction–diffusion models to experimental data of Dictyostelium amoebae reorienting in response to alternating gradients of mechanical shear flow. The experimental readouts we use to fit are spatio‐temporal distributions of a fluorescent reporter for cortical F‐actin labeling the cell front. Experiments performed under different conditions are fitted simultaneously to challenge the models with different types of cellular dynamics. Although the model proposed by Otsuji is unable to provide a satisfactory fit, those suggested by Meinhardt and Levchenko fit equally well. Further, we show that reduction of the three‐variable Meinhardt model to a two‐variable model also provides an excellent fit, but has the advantage of all parameters being uniquely identifiable. Our work demonstrates that model selection and identifiability analysis, commonly applied to temporal dynamics problems in systems biology, can be a powerful tool when extended to spatio‐temporal imaging data. © 2014 The Authors. Published by Wiley Periodicals, Inc.

DIRECTED cell motion is based on three functional modules (i) the formation of cellular protrusions driven by polymerization of actin, (ii) a mechanism to sense extracellular signals, for example, a gradient of chemoattractant, and direct protrusions to the cell front, and (iii) polarization, which is the establishment of a frontrear axis, whereby myosin-II mediates retraction of the cell rear (1)(2)(3). The modular design of cell motility has resulted in it becoming a paradigm of systems biology. In particular, how these modules are integrated to allow cells to navigate in rapidly changing environments has become a focus of theoretical and computational research.
Most models employ a Turing-like (4) local-excitation global-inhibition mechanism, whereby the stronger stimulation of the up-gradient cell end results in local autocatalytic activation of the cell front. At the same time, a fast propagating inhibitory mechanism renders the cell rear unresponsive to stimulation. The theory of reaction-diffusion models is well established and Meinhardt first implemented a model for cell reorientation on a circular domain to study how cells could regain sensitivity at the rear and thus are able to respond to changes in direction of a signaling gradient (5). Most recently, several groups have coupled the Meinhardt model with biophysical models of deformable contours to simulate the deformation and movement of cells in response to a signal gradient (6)(7)(8). Other models have been proposed to address specific questions of signal amplification, sensitivity, and adaptation (9)(10)(11)(12)(13).
Here, we want to compare three representative models for cell reorientation, each of which employs a different regulatory mechanism. The original model by Meinhardt (5), a model by Levchenko and Iglesias (13), and one by Otsuji et al. (12). We are more focused on the problem of model selection, so for their internal workings and their derivation, we refer the reader to an excellent review by Jilkine and Edlestein-Keshet which discusses, among others, all three models in detail (11).
In brief, the Levchenko model has been engineered to achieve perfect adaptation to spatially uniform stimuli, which result in a transient response only, before a new steady-state is achieved. In a gradient, persistent stimulation of a cell front is possible, without the need to temporarily break-down the pattern as in the Meinhardt model. The rationale behind the model by Otsuji et al. (12) is that many signaling components involved in gradient sensing, such as small GTPases of the Rho family are known to exist in an either active or inactive state. Whereas in the models by Meinhardt and Levchenko an increase in signal always causes a stronger response, mass conservation in the Otsuji model takes into account that the total amount of signaling molecules is limited. Specific features of this model are the formation of a strong unique axis of cell polarization and an increased sensitivity at the cell front. The motivation for selecting these models was simply to investigate how more modern models compare to their ancestor, the Meinhardt model. As we wanted to use the same fitting approach for each model we limited ourselves to continuum reaction-diffusion equations, all of which, however, display quite distinct features in their behavior.
In all three models, some of the regulatory mechanisms can be loosely mapped to known biochemical signaling pathways, but all employ a minimal set of regulatory feedback loops, and therefore have a comparatively small number of parameters. This is an important requirement in terms of quantitative modeling that prevents over-fitting and enables selection of structurally identifiable models with unique solutions (14).
Here, we build on our previous work on quantifying actin dynamics in the cortex of moving cells using active contour based methods for cell segmentation and tracking (15)(16)(17). Using fluorescent reporters for polymerized actin as a proxy for cell front activation, we ask: (i) Can we validate different reaction-diffusion models by directly fitting models to time series image data of moving cells? (ii) Will we be able to identify unique sets of parameters?

General Laboratory Reagents
HL5 growth media containing 75 mM glucose was obtained from ForMedium TM (Hunstanton, UK). All other general reagents were purchased from Sigma-Aldrich (St Louis, MO) unless stated.
Experimental Data. The experimental data are fluorescence distributions of a reporter for F-actin (LimED-GFP) in the cortex of Dictyostelium (JH10) cells reorienting in alternating gradients of shear flow as described in (18). Previously, we have shown that the response to shear stress is very similar to that toward a chemoattractant with cells producing a front against the flow. Cells were segmented and tracked using QuimP software [http://go.warwick.ac.uk] (15,19) and fluorescence sampled at 20 equidistant points along the cell cortex. All fluorescence data presented are normalized by dividing through the mean fluorescence in the cell body to account for differences in expression levels, fluctuations in laser intensity and bleaching. Details on microscopy are described in Dalous et al. (18).
Random Motility Experiments. Wild-type Dictyostelium (AX2) cells expressing LimED-RFP were cultured at room temperature in HL5 media containing 75 mM glucose with appropriate antibiotics. Cells were washed twice with KK2 buffer and transferred to glass-bottomed imaging culture plates (Fisher Scientific UK, Loughborough, UK). Actin was visualized using a Personal DeltaVision microscope (Applied Precision, Issaquah, WA) comprising an Olympus UPlanSApo 1003, NA 1.4, oil immersion objective and a Photometric CoolSNAP HQ camera (Roper Scientific, Martinsried, Germany). Captured images were processed by iterative constrained deconvolution using SoftWoRx (Applied Precession) and analyzed using ImageJ (20).

Long Duration Flow Experiments
Wild-type Dictyostelium (AX2) cells expressing ABP120-GFP as a marker for F-actin were cultured in HL5 media containing 75 mM glucose. Cells were washed with KK2 and after 1 h in shaking culture seeded into flow chambers (L 3 W 3 H 50 3 5 3 0.2 mm 3 ), with flow of buffer driven by an air pressure pump system (IB-10902, Ibidi, Martinsried, Germany). A 1 Pa shear flow was applied for 600 s, followed by a 120 s period of no flow. This cycle was repeated once. ABP120-GFP was imaged using a True Confocal Scanner Leica TCS SP5 microscope (Leica Microsystems, Milton Keynes, UK).

Model Fitting
Model fitting and analysis were performed in Potters-Wheel (21), a Matlab toolbox that offers advanced tools for identifiability analyses. The spatial reaction-diffusion models under investigation are systems of coupled, partial differential equations (PDE) with two or three variables. Using a finite difference approximation of the diffusion operator @ 2 C=@x 2 % ðC i21 22C i 1C i11 Þ=ðDxÞ 2 , where C i , C i21 and C i11 denote fluorescence intensities at point i and to the left and right of it, with 1 i N , and N 520 being the number of grid points, Dx the physical spacing between them, the system of PDEs can be transformed into a system of N 3V coupled ordinary differential equations (ODEs), V being the number of variables. The problem is solved on a circular 1D domain with periodic boundary conditions by letting C 0 5C N and C N 11 5C 1 .
When exposed to shear flow cells are aligned with the flow and do not exhibit much variation in shape. We therefore use, for simplicity, an equidistant spacing of fluorescence samples so that Dx is constant. In principle, our approach could be adapted to sampling at arbitrary positions, then requiring interpolation between intensity values at previous time points.
We simulate a spatial profile of shear flow input (source term) using s5ðdycosð2pðn21Þ=N ÞÞ where dy is initially treated as a free parameter which describes the strength of the signal. Later, dy is replaced by an explicit polynomial function dyðPÞ in terms of the absolute shear stress P. Node i5n denotes the up-gradient position on the cell cortex. In the absence of a signal dy is set to zero.
Parameters are fitted by nonlinear least-square minimization using PottersWheel's built-in Trust Region optimization, fitting in logarithmic parameter space (22,23). To avoid getting caught in local minima, models were fitted 10-25 times from different starting points. Where indicated, we fitted to different experimental conditions simultaneously, obtaining a combined v 2 value.

Parameter Profile Likelihood Estimation
To determine the sensitivity and identifiability of parameters, we used PottersWheel's parameter profile likelihood estimation (PLE) tool (24). By iteratively varying a parameter p i about its optimum value and refitting the remaining parameters a v 2 profile of the likelihood for p i can be generated. Where this profile crosses a threshold, v 2 ða; df Þ, lower or upper limits for the confidence interval of p i can be found, here at level a568%. If both limits exist, the parameter is considered identifiable. A value of df 51 yields point-wise confidence intervals, whereas for simultaneous confidence intervals, df equals the number of parameters. The latter is used to determine identifiability.

RESULTS
Our goal is to validate commonly used reaction-diffusion models for cell reorientation by fitting them to time-series image data of cells under well-controlled experimental conditions. The reaction networks of the different models under investigation are depicted in Figure 1A with S denoting the extracellular signal (experimentally we use a gradient of shear flow (18)). Meinhardt's model is based on one autocatalytic activator, A, that produces two inhibitors, B and C (5). The second inhibitor, C, was proposed as an extension to a twovariable model, consisting of only A and B, in order to achieve permanent sensitivity. In the Levchenko model, S promotes simultaneous production of activator A and inhibitor I, which act on a response element R (13). The Otsuji model considers mass conservation of a signaling component, which can be either in an activated form, U, or inactivated, V, whereby S promotes formation of U (12).
Ultimately, activation of a cell front in response to an extracellular stimulus results in formation of a F-actin rich protrusion, which is why we consider F-actin as a reliable readout for front activation. The Dictyostelium cell shown in Figure 1B responds toward shear flow of 2.1 Pa, with the Factin label clearly marking the front facing the flow. Shear flow as signal input can be easily reversed, and the cell can be seen to reorient after changing the flow direction at t 5 0 s. The old front is rapidly degraded, while the new front facing the flow appears at $40 s after flow reversal.
Space-time plots of the cortical fluorescence sampled at 20 nodes allow capturing the entire dynamics in one plot (Fig.  1C). Normalizing the data with respect to the cell circumference has the advantage that data of multiple cells can be averaged, after synchronizing sequences with respect to the time of flow reversal. As illustrated in Figure 1D the population mean of N 5 14 cells responding to 18 flow reversals, provides a clearer indication of F-actin disassembly and reassembly when compared to noisy single cell data. Later, we show that it is possible to fit single cell data, but for our initial model fitting we proceed with population averages. Fitting to single cell data can become computationally demanding when fitting to many cells simultaneously. Figure 1E gives an example of how the fluorescence profile of the averaged cell data along the normalized cell outline is fitted for selected time points using Meinhardt's model. The fitting procedure starts with the flow reversal at t 5 0 s (For details on fitting see Materials and Methods).

The Meinhardt and Levchenko Model Both Fit the Shear-Flow Reorientation Data Well
Having demonstrated how we extract experimental data in form of spatio-temporal maps of F-actin fluorescence and fit them to dynamical models we proceed to compare how each of the three different models could fit to three experimental conditions: reorientation of cells in response to high shear stress of 2.1 Pa; to low shear stress (0.9 Pa); and flow/ no-flow experiments where cells were first oriented under 2.1 Pa, with the flow subsequently being switched off resulting in slow depolarization of cells. These values were chosen since below 0.9 Pa cells respond by making U-turns instead of reversing their orientation and above 2.1 Pa they have difficulties remaining attached to the substratum. We fitted the data using the entire spatial fluorescence profile as shown in Figure  1E, but we summarize the data (Fig. 2) by only plotting the mean fluorescence of the two cell halves (up-gradient and down-gradient). Parameter values and initial conditions are given in Supporting Information Tables S1 and S2. It is apparent that disassembly of F-actin at the old front follows a simple exponential decay under 2.1 and 0.9 Pa with half-lives of T 1/2 5 38 s and 59 s, respectively; assembly of F-actin at the new front is faster under low shear stress of 0.9 Pa where it plateaus after 60 s (Fig. 2). Under the shear stress of 2.1 Pa there is a marked delay of about 30 s before actin polymerization begins leveling around 3 min. The flow/no-flow experiments show slow loss of orientation on the timescale of minutes, but initially there is a slight increase in F-actin after removal of the stimulus. Thus shear flow clearly promotes Factin assembly, but interestingly, at the same time, a higher shear stress slows F-actin assembly. This could be either due to increased mechanical load on the F-actin network at higher shear stresses, or negative feedback in the biochemical signal transduction pathway.
Shear flow affects the time-scales of loss and gain of fluorescence in intricate ways. Initially, we avoided dealing with absolute values of shear stresses as signal input, and let the

Identifiability Analysis
The model fitting described previously provides parameters that maximize the goodness of fit, however, it is important to consider if these are the only combination of parameters that can explain the data. To address this, we have used profile likelihood estimations, where values for one particular parameter are changed over a defined region and the model is repeatedly refitted to compute how changing the parameter affects the goodness of fit. Identifiable parameters are characterized by a parabolic v 2 profile which indicates a unique, optimal parameter value (Fig. 3A). Unidentifiable parameters do not affect the quality of the fit, and consequently have a flat profile; they might however be constrained by an upper or a lower limit. Unidentifiability can be linked either to the structure of a model, requiring changes to the model itself or a lack of quality experimental data causing practical unidentifiabilities. The profile likelihoods for the Meinhardt model show that all but one parameter were identifiable with r b practically non identifiable. In the Levchenko model four out of ten parameters were practically unidentifiable.

A Simplified Two-Variable Version of the Meinhardt Model is Uniquely Identifiable
As shown previously, both, the models by Meinhardt and Levchenko were not fully identifiable. Different approaches exist to make models identifiable, changing either the model or the experimental design (24). In the Meinhardt model, we observed that the first inhibitor, B, remained close to 1, thus we first considered reducing the model to two variables, setting @B=@t to zero. Ideally, it would still capture the same dynamics of A but become identifiable as the previously unidentifiable parameter r b was dropped. However, we were unable to find a single value for B around 1, nor a simple linear expression in terms of the external shear stress P, which fitted all three experimental conditions simultaneously. We, therefore, determined optimum values of B for each condition and fitted a quadratic, BðPÞ511 b 0 ðP 2 1b 1 PÞ, in terms of P, which made it possible to simplify dy to dy(0) 5 0, and dyðPÞ5const50:0128. Contrary to the original three-variable Meinhardt model, the reduced model has the advantage that it depends explicitly on the external shear stress. In principle, the exact nature of B(P) could be tested through fitting additional experiments with different shear stresses P, which however is beyond the scope of the current study.
We next performed profile likelihood estimations (Fig.  4A) for multi-experiment fits of the two-variable model to the same three experimental conditions as in Figure 2. The twovariable model generated a similar good fit (v 2 5184 compared to v 2 5178 for the three-variable model, half-lives for the loss of F-actin at the old front are almost identical). Moreover the parameter sets obtained when fitting were also similar between the two and three-variable model (Supporting Information Table S1 and Fig. 4B).

Using the Models to Make Predictions About the Persistence of Front Activation
Both, the Meinhardt and the Levchenko models fitted shear flow reversals reasonably well. We wanted to test whether we would be able to predict the outcome of a new experiment Figure 2. Actin relocalization data enables comparison of three proposed models of cell polarity: to provide a comparison of repolarization cells were split into two halves such that the front (solid line) was determined using the average cortical fluorescence obtained from the nodes 6-15 as illustrated in Figure 1B. The average cortical fluorescence for the back half of the cell (dotted line) was determined from nodes 1-5 and 16-20. Biological data was obtained as described in (18) for a high shear stress (P 5 2.1 Pa), a low shear stress (P 5 0.9 Pa) and cessation of flow. Mean data for the high shear stress was obtained from 18 responses from 14 cells, for the low shear stress 10 responses from 5 cells was analyzed, while 13 responses from 9 cells are shown in the cessation data set. For all data the shaded area represents standard error of the mean. Models as defined in Figure 1 were simultaneously fitted to all biological data sets using Pot-tersWheel (see Methods). All parameters were conserved between data sets with the exception of dy, the asymmetry of the external stimulus. The goodness of fit, v 2 , is 178, 155, and 361 for the Meinhardt, Levchenko, and Otsuji models, respectively.
Original Article using the sets of parameters obtained previously. We exposed a cell to two cycles of flow-induced (1 Pa) polarization and observed a persistent migration of the cell toward the flow source (Fig. 5A). Cessation of the flow for 2 min caused a loss of polarity. We performed simulations where we first wanted to qualitatively reproduce the persistent activation of the cell front as seen by the steady migration of the example cell for 10 min (Fig. 5B). These were initialized with uniform conditions and the system equilibrated, before replicating the signal behavior detailed above. Using the parameter set obtained in Figure 3, the Levchenko model failed to produce a stable front. To obtain a stable front a rather drastic change to the model was required, for example, increasing D 1 by four orders of magnitude and k 2A by a factor of 10. (Supporting Information Table S1 and Fig. 5C). Using the parameter set obtained in Figure 4 for the modified Meinhardt model, a single front was obtained but it rapidly split into three. A single persistent front could be obtained by reducing diffusion of the inhibitor C by 20% (Fig.

Fitting Spontaneous Front Activation in Randomly Migrating Cells
Previously, we aggregated and synchronized data from multiple cells to obtain averaged data where noise is significantly reduced. Next, we attempted to simultaneously fit single cell data of randomly migrating Dictyostelium cells to determine how the models performed on more complex and noisy data (Figs. 6A and 6B). We fitted our biological data using both the Meinhardt and Levchenko models. Whereas the models are deterministic, the observed random patterns are clearly driven by noise. By treating the initial activator concentrations as free parameters, we only account for noise at the start of the time series.
We observed that for the Meinhardt model four parameters were an order of magnitude lower when compared to parameters obtained in Figure 5 (Supporting Information  Table S1), three were of the same order and D C could essentially be set to zero. Parameters in the Levchenko model also vary greatly between reorientation and random motility experiments (Supporting Information Table S1).
Both models captured some of the intrinsic dynamics surprisingly well (Figs. 6A, and 6B): (i) a front which abruptly disappears, (ii) a front which splits into two, and (iii) a broader low intensity F-actin crescent at the cell rear. Regions IV and VI of a second cell (Fig. 6B) again resemble front splitting, while region V denotes the merger of two fronts. Fitting reaction-diffusion models to image data has been successfully applied in image enhancement in many areas (25), Our example shows that the same concepts, usually applied to single images, can be extended in a straight-forward manner to filtering time series image data, thus aiding the model-based analysis of complex stochastic time-series data.

DISCUSSION
A number of sophisticated mathematical models that couple models for cell orientation to cell deformation have recently been published (6)(7)(8). Interestingly, they all employ the original model for cell orientation by Meinhardt (5), which as we have shown demonstrates good agreement with experimental data of Dictyostelium cells responding to three different experimental conditions (18). Significantly, the model is able to make reasonable predictions of cortical F-actin dynamics during cell reorientation for up to 2 min, which is remarkable given that the turnover of the entire F-actin system in Dictyostelium is on the timescale of seconds (26,27).
Chemotactic receptors have been very well characterized but only recently has light been shed on putative mechanosensors, in particular PKD2 Ca 21 channels (28). Previously we Original Article have shown that in Dictyostelium cells reversal of cell polarity in response to a shear flow reversal is very similar to the response seen when reversing a chemotactic gradient. This suggests that both sensor systems might be linked to one common pathway regulating actin polymerization and cell polarization (18). It is interesting to note that the models originally developed to investigate chemotaxis are equally applicable to mechanotaxis. The estimation of model parameters from image data and identifiability analysis in the context of diffusion processes is an emerging area of research (29). We have reduced the Meinhardt model to a two-variable model, which is uniquely identifiable. Although there is no direct correspondence with known biochemical pathways, the activator variable in the model captures, remarkably well, the dynamics of F-actin assembly. The parameters we determined yield biochemically realistic timescales of cell front activation and repolarization, which will be appreciated by modelers trying to build more complex models that integrate actin dynamics and protrusive behavior.
The behavior of the reduced two-variable Meinhardt model is almost identical to the three-variable model. Although the second inhibitor in the Meinhardt model is often regarded as an improved extension of a model with only one inhibitor, adding permanent sensitivity, its local action can completely replace the first global inhibitor. In line with Figure 5. Simulated responses to long duration stimulation: (A) image sequence of wild-type (AX2) cells expressing ABP120-GFP (green). Cells were imaged for a total period of 1,560 s. Following an initial rest period (120 s) cells were exposed to a continuous flow for a duration of 600 s. The cycle was then repeated for a further 720 s. Arrows indicate presence and the direction of a shear stress P 5 1 Pa; (B) Position of the cell front (relative to the bottom of the image) throughout the time course. Asterisks denote times of the five individual images shown. Scale bar, 10 lm; (C) simulations of the biological data in (A) using the Levechenko model. The arrow denotes the position where the external signal is the strongest. Adjustment of model parameters D I and k -A in the Levchenko model allows a persistent front to be obtained; (D) simulations of the biological data in (A) using the reduced Meinhardt model as described in Figure 4 results in breaking up of a single front; and (E) decreasing the diffusion constant for inhibitor C from 9.768 3 10 22 to 7.064 3 10 22 produces a persistent front. All models were started with uniform conditions and simulated in the absence of signal for 5,000 s before the signal input from experiment was replicated. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.] that we found that the first inhibitor remained constant over time, justifying its removal as an independent variable. Similarly the Levchenko model achieves the same behavior using only one inhibitor, and also captures many aspects of F-actin dynamics during reorientation and random motility.
We have shown that care is needed when trying to make predictions as initially we were unable to achieve long-term persistence of front activation. Both the Meinhardt and Levchenko models required changes in their parameters. Significantly, this helped to constrain the models further; indeed the Meinhardt model was still able to fit our original set of data under the new constraint.
Randomly migrating cells display several competing fronts, which aids the probing of their environment by increasing the sampling frequency. The presence of a strongly orientating stimulus such as a chemotactic agent or shear flow requires cells to enter a different state where only a single front survives. Thus, it is expected that this new state is reflected by changes in parameters. In theory, going from several fronts to a single one is equivalent to increasing the wavelength of the pattern, which can be achieved by increasing activator diffusion. Here we confirm, for the Meinhardt model, that the diffusion of the activator in randomly migrating cells is by a factor of 20 lower, when compared to the shear flow experiments. The rate of inhibitor Original Article diffusion becomes essentially negligible, suggesting a qualitative change in the model such that there is no lateral inhibition of activated peaks any more. This lack of lateral inhibition results in patterns of activation which are not evenly spaced, which is in fact a notable feature of random migration. Therefore, the changes in fitted diffusion rates are in agreement with the experimental observations. At the same time, fronts in randomly migrating cells are less stable, which explains why some of the kinetic parameters are also required to change. Interpretation of the exact changes is however difficult. Although the rates of production and decay of the local inhibitor C are, for example, reduced by a factor of 4, its local concentration can still increase rapidly, because it does not diffuse. This might contribute to the observed shorter lifetime of activated fronts in randomly migrating cells. Recently, membrane tension has been discussed as an inhibitory mechanism restricting the growth of protruding fronts, which because of its physical nature could explain very fast diffusive spreading (30,31). Therefore, one possibility supported by the change in model parameters could be that cells globally increase membrane tension when switching from random migration to a strongly polarized mode of movement.
Similarly to our example of modeling cell orientation, inverse modeling of spatio-temporal cellular dynamics has been employed in the context of photobleaching or photoactivation experiments to study the mobility, mass transport, or binding of cellular constituents (32)(33)(34). Perturbing the intrinsic dynamics by photobleaching or activation could nicely complement our approach, thereby trying to match diffusion in the model with the mobility of molecular players known to be involved in signaling to the actin cytoskeleton. More detailed biological models exist for receptor/G-protein networks, and signaling to the downstream modules controling polarity/myosin-II contractility and actin reorganization (35). Recently, Skoge et al. (36) have developed a model for Ras activation in chemotaxing Dictyostelium cells which includes memory effects to explain why cells do not reverse direction in the wake of a wave of chemoattractant. The model consists of seven equations with 24 parameters and has been successfully fitted to activated Ras levels measured at the front and the back of cells. We believe that our proposed framework will be a valuable tool to compare this and other recent models in the future. For example, to learn what degree of complexity is required to explain particular experimental findings which cannot be easily explained by the very simple models discussed in the current paper.