Estimating cell diffusivity and cell proliferation rate by interpreting IncuCyte ZOOMâ„¢ assay data using the Fisher-Kolmogorov model
IncuCyte ZOOMâ„¢ Assay
We perform a monolayer scratch assay using the IncuCyte ZOOMâ„¢ live cell imaging system
(Essen BioScience, MI USA). This system measures scratch closure in real time and
automatically calculates the relative wound density within the initially-vacant area
at each time point. The relative wound density is the ratio of the occupied area to
the total area of the initial scratched region. All experiments are performed using
the PC-3 prostate cancer cell line 16], which is obtained from the American Type Culture Collection (ATCC, Manassas, USA).
Cells are routinely propagated in RPMI 1640 medium (Life Technologies, Australia)
in 10 % foetal calf serum (Sigma-Aldrich, Australia), with 110 u/mL penicillin, 100
?g/mL streptomycin (Life Technologies), in plastic flasks (Corning Life Sciences, Asia
Pacific) in 5 % CO
2
and air in a Panasonic incubator (VWR International) at 37 °C. Cells are regularly
screened for Mycoplasma (ATCC). Cells are removed from the monolayer using TrypLEâ„¢(Life Technologies) in
phosphate buffered saline, resuspended in medium and seeded at a density of 20,000
cells per well in 96-well ImageLock plates (Essen BioScience). After seeding, cells
are grown overnight to form a spatially uniform monolayer. We use a WoundMakerâ„¢(Essen
BioScience) to create uniform, reproducible scratches in all the wells of a 96-well
plate. After creating the scratch, the medium is aspirated and the wells are washed
twice with fresh medium to remove any cells from the scratched area. Following the
washes, for the control assay, 100 ?L of fresh medium is added to each well. We also perform a series of experiments where,
following the washes, fresh medium containing different concentrations of EGF (Life
Technologies) is added to the wells. The concentrations of EGF we use are: 25, 50,
75, 100 and 125 ng/mL. We will refer to these assays as EGF-25, EGF-50, EGF-75, EGF-100
and EGF-125, respectively. Once the fresh medium is added, the plate is placed into
the IncuCyte ZOOMâ„¢ apparatus and images of the collective cell spreading are recorded
every 2 hours for a total duration of 46 hours. For the control assay and each different
EGF concentration we perform three identically prepared experimental replicates (n=3).
Image analysis
We use Matlab’s Image Processing Toolbox 23], 24] to estimate the position of the leading edge of the spreading cell population in
the IncuCyte ZOOMâ„¢ images. The experimental image is imported and converted to greyscale
using the imread and rgb2gray commands, respectively. We detect edges in the images using edge with the Canny method 25] and automatically-selected threshold values. Detected edges outside of these threshold
values are ignored. Remaining edges are dilated using the imdilate command and a structuring element, defined using strel, with a circular element of size 15. Using the bwareaopen command with a component size of 10,000 pixels, we remove any remaining vacant spaces
in the image while preserving the vacant scratch. Edge dilation is reversed using
the imerode command with the same structuring element defined previously to erode the image. Finally,
edges within the image are smoothed using medfilt2 and the area of the remaining vacant space, A(t), representing the vacant area, is estimated using the regionprops command.
We calculate the position of the leading edge, which we define to be the distance
between the centre of the experimental domain and the position of the leading edge
using
(1)
where L x
is the horizontal width of the image and L y
is the vertical height of the image. For all experiments we have L x
=1970?m and L y
=1430?m. Eq. (1) allows us to examine the time evolution of the scratched area in terms of L E
(t), which is the half-width of the scratch (Fig. 1(a)).
Mathematical model
We interpret our experimental results using the Fisher-Kolmogorov equation 26]–28], which is a continuum reaction-diffusion model describing the spatiotemporal evolution
of cell density in a population of cells where cell migration is driven by random
(undirected) cell motility and cell proliferation is driven by carrying capacity limited
logistic growth. The Fisher-Kolmogorov equation, and extensions of the Fisher-Kolmogorov
equation, have been previously applied to in vitro29]–31] and in vivo32], 33] data describing collective cell spreading in a range of contexts including wound
healing 34], 35], tissue repair 3], 4] and malignant spreading 8]–10], 36], 37].
Since our scratch assay takes place on a two-dimensional substrate (Fig. 1(a)–(c)), we start with the two-dimensional analogue of the Fisher-Kolmogorov equation in
Cartesian coordinates
(2)
where [cells/ ? m2
] is the cell density, or average number of cells per unit area, at location (x,y) and time t. For our experiments we have 0?x?L x
and 0?y?L y
. There are three parameters in the Fisher-Kolmogorov equation: (i) the cell diffusivity,
D [ ? m2
/h], (ii) the cell proliferation rate, ? [/h], and (iii) the carrying capacity density, K [cells/ ? m2
]. The proliferation rate, ?, is related to the cell doubling time t d
=log
e
(2)/?. We note that we make the standard assumption that D, ? and K are constants 1]–4], 29], 34]
Since the initial cell monolayer is spatially uniform and the initial scratch is made
perpendicular to the x-direction (Fig. 1(a)), we can simplify the mathematical model by averaging in the y-direction 8]–10]. To do this we average the two-dimensional cell density
(3)
which allows us to write Eq. (2) as a one-dimensional partial differential equation
(4)
In general, approximating a two-dimensional nonlinear partial differential equation,
such as Eq. (2), by an averaged one-dimensional nonlinear partial differential equation, such as
Eq. (4), can introduce an averaging error. However, for initial conditions such as ours
where the initial density is independent of the vertical direction, this error vanishes,
and a detailed analysis of this error is presented elsewhere 38], 39]. The initial condition for Eq. (4) is given by the width of the scratch (Fig. 1(a))
(5)
where C0
is the initial density of cells in the monolayer and L E
(0) is the initial position of the leading edge. Since we use a WoundMakerâ„¢tool to
create uniform scratches in all experimental replicates, the initial condition, given
by Eq. (5), applies to all experiments and cannot be varied.
The physical distribution of cells in each experiment extends well-beyond the L x
×L y
rectangular region imaged by the IncuCyte ZOOMâ„¢ apparatus. Therefore, since the cells
are spatially uniform except for the scratched region, there will be no net flux of
cells across the vertical boundaries along the lines x=0 and x=L x
. We model this by using zero-flux boundary conditions
(6)
These boundary conditions do not imply that cells are stationary at x=0 and x=L x
. Instead, these boundary conditions imply that the cell density profile is spatially
uniform, ? C(x,t)/? x=0, so that there is no net flux of cells across the boundaries at x=0 and x=L x
.
We solve Eq. (4) using a finite difference numerical method 40]. The spatial domain, 0?x?L x
, is discretised uniformly with grid spacing ? x, and the spatial derivatives are approximated using a central-difference approximation
40]. This leads to a system of coupled nonlinear ordinary differential equations that
are integrated through time using a backward-Euler approximation with constant time
steps of duration ? t40]. The resulting system of coupled nonlinear algebraic equations are linearised using
Picard (fixed-point) iteration, with absolute convergence tolerance ?41]. The associated tridiagonal system of linear equations is solved using the Thomas
algorithm 40]. For all results presented here we always chose ? x, ? t and ? so that our numerical algorithm produces grid-independent results.
We also apply Eq. (4) to some simplified situations where we focus on the time evolution of the cell density
in small subregions, located well-behind the initial position of the scratch, where
the cell density is approximately spatially uniform. This implies that C(x,t)?C(t) within these subregions 18], 21], 22]. Since the cell density is approximately spatially uniform we have ? C(x,t)/? x=0, and the first term on the right of Eq. (4) vanishes and, subsequently in these subregions, the partial differential equation
simplifies to the logistic equation,
(7)
whose solution is given by
(8)
where C(0)=C0
is the initial density at t=0. The simplification of approximating Eq. (4) by Eq. (7) in subregions well-behind the leading edge where the cell density is spatially uniform
does not imply that cells are stationary in these subregions. Instead, Eq. (7) represents the situation where there is no gradient in cell density and cells are
free to move amongst the extracellular space within these subregions. The key advantage
of applying this approximation is that cell motion in these spatially uniform subregions
does not contribute to any temporal changes in cell density. Instead, when the cell
density is spatially uniform, any temporal change in cell density is solely associated
with the proliferation term in Eq. (4) 18], 21], 22].
Parameter estimation
We estimate the three parameters in the Fisher-Kolmogorov model using a sequential
approach. First, using cell counting, we estimate the parameters governing cell proliferation:
K and ?. Second, using data describing the temporal changes in the position of the leading
edge, we estimate the cell diffusivity, D. Although it is possible to use a different approach, based on a multivariate regression
technique to estimate D, ? and K simultaneously, we prefer to estimate these parameters sequentially. Estimating the
three parameters sequentially, one at a time, emphasises the differences in the interpretation
of these parameters, as well as emphasising the differences in the mechanisms of cell
proliferation and cell motility. If, instead, a multivariate approach is used to estimate
the three parameters simultaneously, we anticipate that the interpretation of the
mechanisms associated with these parameters might not be obvious as it is in our approach.
Carrying capacity density
To estimate K we focus on experimental images from the latter part of the experiment, t=46 h, where the cell population has grown to confluence. We identify three smaller
subregions, located well-behind the initial leading edge, and count the number of
cells within each subregion, N. To quantify the variability in our estimate we analyse three different subregions
in each image and count N in each replicate subregion. Using this data we estimate the average carrying capacity
density as
(9)
where ?N? is the average number of cells within the subregion of area A SR
=3.789×10
4? m2
. To examine whether EGF has any impact on the carrying capacity density we estimate
K for the control assay and for each experiment treated with a different EGF concentration.
Figure 2 shows IncuCyte ZOOM™ images at t=46 h with the location of three subregions superimposed. The images in Fig. 2(a)–(c) show the control, EGF-50 and EGF-100 assays, respectively. We note that the location
of all three subregions in each image is located well-behind the initial position
of the scratch (Fig. 1(a)) so that after t=46 h the local density of cells within each subregion has grown to confluence. To
quantify the variability in our estimate of K, we calculate the sample standard deviation for each EGF concentration, and report
results as a mean value for K, with the variation in our estimate given by plus or minus one standard deviation
about the mean. Results are summarised in Table 1.
Fig. 2. Final time experimental images (t=46 hours) for three IncuCyte ZOOMâ„¢ assays for (a) Control, (b) EGF-50, and (c) EGF-100. The three coloured boxes indicate the location of the three subregions
used to estimate K and ?. Each coloured square within the subregions indicates the centre of an individual
cell in the cell counting step. Scale bar corresponds to 300 ?m
Table 1. Estimated K, ?, D and C0
values for PC-3 cells for different EGF concentrations. Results are reported as a
mean, with the estimate of the variability given in the parenthesis
Proliferation rate
The logistic equation, given by Eq. (7), describes the time evolution of the cell density where there is, on average, no
spatial variation in cell density. To apply the logistic equation to our data we analyse
three subregions within each IncuCyte ZOOMâ„¢ image at several time points during the
assay. Counting the total number of cells in each subregion and dividing by the area
of the subregion gives an estimate of the local cell density in that subregion. In
all cases the subregions we considered always started off with approximately 20–30
cells at t=0 h. Repeating this procedure for three different subregions, at fixed locations,
for each experimental replicate, at five different time points, allows us to calculate
the average cell density as a function of time, C(t). With this data, together with our previous estimates of K, we find the value of ? in Eq. (8) that matches our C(t) data across several time points. For consistency, when we estimate ? we always analyse the same three subregions that we used previously to estimate K. The location of these three subregions is shown in Fig. 3. We estimate the initial cell density, C(0), from the first image taken immediately after the scratch is made at t=0 h. Images of the assay in Fig. 3(a)–(e) correspond to 0, 8, 16, 24 and 46 h, respectively. The location of each subregion
is chosen to be well-behind the initial position of the leading edge of the population
so that the cell density is approximately spatially uniform locally within each subregion.
In each subregion C(t) increases with time, and we attribute this increase to cell proliferation. The data
in Fig. 3(f) shows the time evolution of the average cell density, C(t), calculated by averaging the three estimates of cell density from each subregion,
at each time point. Using our previous estimate of K, we estimate ? by matching the solution of Eq. (7) with the observed C(t) data.
Fig. 3. a–e Time evolution of an EGF-75 IncuCyte ZOOM™ assay. Images taken after (a) 0, (b) 8, (c) 16, (d) 24, and (e) 46 h after the scratch was performed. The three coloured boxes indicate the location
of the three subregions used to calculate K and ?. Each coloured square within the subregions indicates the centre of an individual
cell in the cell counting step. Scale bar corresponds to 300 ?m. f Comparison of the average experimental cell density C(t) (crosses) and the logistic growth curve using our estimates of K and ? (solid)
For each EGF concentration we have three sets of data describing the temporal variation
in average cell density per experimental replicate. For each set of time series data
we use Matlab’s lsqnonlin function, a nonlinear least squares minimisation routine 42], to estimate ?. To quantify the average proliferation rate we calculate the average ? value by averaging the three estimates from each experimental replicate. A comparison
of the resulting logistic growth curve using our average estimate of K and ? with the observed C(t) data is given in Fig. 3(f) for the EGF-25 experiment, indicating that the solution of Eq. (7) matches the data reasonably well. To quantify the variability in ? we report our mean estimate of ? and the variation as the mean plus or minus one standard deviation. We also report
the mean and variability in the C0
values used to obtain estimates of ?. Results are summarised in Table 1.
Cell diffusivity
Several different approaches have been used in previous studies to estimate D from in vitro assays describing collective cell spreading processes. For example, Treloar et al.
22] estimate D in a circular barrier assay by applying two different methods to the same data set.
First, they estimate D using a cell labelling and cell counting technique to provide an estimate of the
cell density profiles near the leading edge of the spreading population. Treloar et
al. 22] calibrate the solution of a partial differential equation to that data to give an
estimate of D which matches the position and shape of the spreading cell density profile. Second,
using the same data set, Treloar et al. 22] use automated leading edge image analysis 23], 24] to quantify temporal changes in the position of the leading edge of the spreading
cell density profile without counting individual cells. Treloar et al. 22] calibrate their model to this leading edge data to obtain a second estimate of D. Given that the two approaches implemented by Treloar et al. 22] produce similar estimates of D, here we choose to estimate D using a similar technique based on leading edge data since this is the most straightforward
approach which avoids the need for labelling and counting individual cells within
the spreading population.
Figures 4(a)-(d) show IncuCyte ZOOMâ„¢ images at 0, 10, 20 and 30 h, respectively. The position of
the two detected leading edges of the spreading population is superimposed on each
image. A visual comparison of how the position of the detected leading edges changes
with time suggests that the initially-vacant region closes symmetrically with time.
The edge detection results allow us to calculate the area of the vacant region, A(t), and with this information Eq. (1) allows us to estimate the half-width, L E
(t), which is decreasing function of time. Previously, Treloar et al. 22], 24] showed that the location of the automatically detected leading edge corresponds to
a cell density of approximately 2 % of the carrying capacity density.
Fig. 4. a–d Indicate the area of remaining vacant space, A(t), as determined by the edge detection algorithm at a 0, b 10, c 20, and d 30 h for the control assay. The position of the detected leading edge is given in
green. The straight vertical lines superimposed on a (white) indicate the average width of the scratch, 2L E
(t). Scale bar corresponds to 300 ?m. e Average L E
(t) data estimated from the control assay experimental images (blue). The error bars correspond to one standard deviation about the mean. Numerical L E
(t) data (red), corresponding to the numerical solution of Eq. (4) using the relevant estimates of D, ? and K (Table 1). (f) Evolution of C(x,t) profiles at t=0, 10, 20, 30 h corresponding to the numerical solution of Eq. (4) using the relevant estimates of D, ? and K (Table 1). Arrows indicate the direction of increasing time. Numerical solutions of Eq. (4) correspond to ? x=1?m, ? t=0.1 h and ?=1×10
?6
. The vertical lines show the locations of the subregions where the estimates of ? and K were obtained
Given our previous estimates of K and ?, and assuming that the position of the detected leading edge corresponds to the location
where the density is 2 % of the carrying capacity, we use Matlab’s lsqnonlin function to find an estimate of D that minimises the difference between the observed time series of L E
(t) and the time series of L E
(t) data from the numerical solution of Eq. (4). We present an example of the match between the experimental measurements of L E
(t) and numerical prediction of L E
(t) in Fig. 4(c). For all time points, the numerical estimate of L E
(t) is always within one standard deviation of the mean of the experimental measurements.
Given our estimates of D, ? and K, we can use our numerical solution of Eq. (4) to explore how C(x,t) varies across the entire width of the domain, for the duration of the assay, as
illustrated in Fig. 4(f). These profiles show that the cell density remains approximately spatially uniform
well-behind the initial location of the scratch. In fact, we have indicated the position
of the location of the various subregions used to estimate K and ? on the profiles in Fig. 4(f), and we see that the predicted cell density profile is spatially uniform, ? C(x,t)/? x=0, at these locations for the duration of the assay, which is visually consistent
with the subregions presented in Fig. 2. The cell density in the three subregions well-behind the initial location of the
scratch increases with time owing to cell proliferation. These profiles also show
the cell density front near the location of the scratch moves inward to close the
initially-scratched region with time.
Using our approach we calculate an average value of D by estimating the cell diffusivity for each experimental replicate and then averaging
the results. We note that Treloar et al. found that varying the Matlab edge detection
threshold parameters led to a small variation in the position of the detected leading
edge corresponding to a cell density in the range of approximately 1–5 % of the carrying
capacity density 24]. To quantify the variability in our estimate of D we repeated the edge detection by assuming that the position of the detected leading
edge corresponds to both 1 and 5 % of the carrying capacity density. Results are summarised
in Table 1. We note that the maximum value of D is obtained by assuming that the detected leading edge corresponds to 5 % of the
carrying capacity since this upper bound implies additional spreading. Conversely,
the minimum value of D is obtained by assuming that the detected leading edge corresponds to 1 % of the
carrying capacity since this lower bound implies less spreading.