In the above, we saw that a Brownian motion model underestimated the
MSD of the real T-cell dataset. Here, we show how we can do better by
fitting models on the MSD.
5.1 Before we start: dimensionality and imaging windows
Before we start, let’s look at the data we are fitting in more
detail. If we look at the dimensions of the T-cell dataset:
## t x y
## min 24 0.4 0.400
## max 960 409.2 334.724
we see that the imaging window is limited to roughly 200 x 200 x 50
\(\mu\)m; thus, it is smaller in the
z-dimension.
A limited imaging window like this can cause artefacts in the MSD
curve because cells can not be followed after they leave the imaging
window. Thus, for larger \(\Delta t\),
we underestimate the mean (squared) displacement because the faster
cells have left the imaging window and do not end up in our data. In our
case, this will be particularly problematic in the z-dimension. Thus,
for now we will not fit the 3D data, but instead use the 2D projection
on the xy-plane:
Txy <- projectDimensions( TCells, c("x","y") )
Let’s look at the MSD curve we get:
# Compute MSD
msdData <- aggregate( Txy, squareDisplacement, FUN = "mean.se" )
# Scale time axis: by default, this is in number of steps; convert to minutes.
tau <- timeStep( Txy ) / 60 # in minutes
msdData$dt <- msdData$i * tau
# plot
ggplot( msdData, aes( x = dt, y = mean ) ) +
geom_ribbon( alpha = 0.3, aes( ymin = lower, ymax = upper ) ) +
geom_line() +
labs( x = expression( Delta*"t (min)" ),
y = expression( "displacement"^2*"("*mu*"m"^2*")") ) +
coord_cartesian( xlim = c(0,NA), ylim = c(0,NA), expand = FALSE ) +
theme_bw()
Two things are worth noting here:
the larger \(\Delta\)t, the
higher the uncertainty (SE); this is because at larger \(\Delta\)t, we have fewer subtracks of that
length, and thus less data to compute a mean on;
the MSD curve has a weird shape. This is actually the case
because of the limited imaging window problem described above. We will
therefore use only the first part of the curve (the first 5 minutes),
where we have enough data and the imaging window problem remains
manageable.
fitThreshold <- 5 # minutes
Now, let’s see how we can fit the two simulation models:
brownianTrack and beaucheminTrack on this MSD
curve.
5.2 Fitting Brownian motion based on the diffusion coefficient
The brownianTrack function takes 4 arguments: the number
of steps to simulate, the number of dimensions to simulate in, and the
mean and sd of the step sizes $(x, y, z) in each dimension. The number
of steps depends on the time we wish to simulate, the dimension is 2
(since we projected the T cells on the xy plane) and the mean of an
unbiased random walk is always 0. That leaves the standard
deviation.
For Brownian motion, the dimension-wise variance of the step size is
related to the “diffusion coefficient” \(D\) as:
\[\text{var} = 2D\tau\] with \(\tau\) the step duration. For the standard
deviation, we get: \[\sigma_\text{dim} =
\sqrt{2D\tau}\] we already know the \(\tau\) of our data, and can obtain \(D\) by fitting Furth’s formula for the mean
squared displacement \(\langle d^2
\rangle\):
\[\langle d^2 \rangle (\Delta t) = 2Dn_d
\big( \Delta t - P(1-e^{-\Delta t/P} ) )\] Here, \(D\) is the diffusion coefficient, \(P\) the persistence time, \(\Delta t\) the time resolution we are
looking at, and \(n_d\) the number of
dimensions.
We can fit this function using R’s nls() on the MSD data
to get \(D\) and \(P\):
fuerthMSD <- function( dt, D, P, dim ){
return( 2*dim*D*( dt - P*(1-exp(-dt/P) ) ) )
}
# Fit this function using nls. We fit only on the data where
# dt < fitThreshold (see above), and need to provide reasonable starting
# values or the fitting algorithm will not work properly.
model <- nls( mean ~ fuerthMSD( dt, D, P, dim = 2 ),
data = msdData[ msdData$dt < fitThreshold, ],
start = list( D = 10, P = 0.5 ),
lower = list( D = 0.001, P = 0.001 ),
algorithm = "port"
)
D <- coefficients(model)[["D"]] # this is now in units of um^2/min
P <- coefficients(model)[["P"]] # persistence time in minutes
D
## [1] 19.91724
## [1] 0.3069776
Note that we will not use \(P\).
However, we fit Furth’s formula with a non-zero persistence time \(P\) because we would otherwise
underestimate \(D\).
Now, given \(D\), simulate some
tracks:
# simulate tracks using the estimated D. We simulate with time unit seconds
# to match the data, so divide D by 60 to go from um^2/min to um^2/sec.
# The dimension-wise variance of brownian motion is 2*D*tau (with tau the step
# duration), so use this to parametrise the model:
tau <- timeStep( Txy )
brownianScaled <- function( nsteps, D, tau, dim = 2 ){
tr <- brownianTrack( nsteps = nsteps, dim = 2, sd = sqrt( 2*(D/60)*tau ) )
tr[,"t"] <- tr[,"t"]*tau # scale time
return(tr)
}
brownian.tracks <- simulateTracks( 1000,
brownianScaled( nsteps = maxTrackLength(Txy ),
D = D, tau = tau ) )
plot( brownian.tracks[1:10], main = "brownian motion")
Finally, get MSD and compare it to data:
# get msd of the simulated tracks and scale time to minutes
msdBrownian <- aggregate( brownian.tracks, squareDisplacement, FUN = "mean.se" )
tau <- timeStep( brownian.tracks ) / 60 # in minutes
msdBrownian$dt <- msdBrownian$i * tau
# Get Furth MSD using fitted P and D for comparison
msdFuerth <- data.frame( dt = msdBrownian$dt,
mean = fuerthMSD( msdBrownian$dt, D, P, dim = 2 ))
# plot
ggplot( msdBrownian, aes( x = dt, y = mean) ) +
geom_point( data = msdData, shape = 21, aes( fill = dt < fitThreshold ), color = "blue", show.legend = FALSE ) +
geom_line( data = msdFuerth, color = "red", lty = 2 ) +
geom_line( )+
scale_fill_manual( values = c( "TRUE" = "blue", "FALSE" = "transparent" ) ) +
scale_x_log10( expand = c(0,0) ) +
scale_y_log10( expand = c(0,0) ) +
labs( color = NULL, x = expression( Delta*"t (min)"),
y = expression( "displacement"^2*"("*mu*"m"^2*")") , fill = NULL ) +
theme_bw()
Here, we see the data in blue, where the colored circles are the data
the fit is based on. The red line is the Furth’s MSD, the black line is
the Brownian model. The Brownian model does not describe the data well
especially on short time scales, because it does not contain the
persistence that the Furth formula does. Nevertheless, this fit is
better than the one with matched displacements on short time scales.
Next, let’s look at a slightly more complicated model.
5.3 Fitting the Beauchemin model
We now use a similar approach to fit the model proposed by Beauchemin et al
(2007), which was designed for T-cell migration in the lymph
node.
The MSD of this model can be derived analytically (Textor et al (2013),
proposition 12) :
\[\langle d^2 \rangle (\Delta t)= 2Mt -
2Mt_\text{free} \cdot \begin{cases}
\frac{1}{3} & t \geq t_\text{free} \\
\frac{1}{3} \big(t/t_\text{free}\big)^3 - \big(t/t_\text{free}\big)^2 +
\big(t/t_\text{free}\big) & t < t_\text{free}
\end{cases}\]
where \(M\) is the motility
coefficient, which should be the same as our diffusion coefficient \(D\) (see above). \(M\) is related to the model parameters (Textor et al (2013),
equation 2):
\[M = \frac{v_\text{free}^2\cdot
t_\text{free}^2}{6(t_\text{free}+t_\text{pause})}\]
We should note here that the Beauchemin has three parameters (\(v_\text{free}\), \(t_\text{free}\), \(t_\text{pause}\)), but that these are
underdefined based on MSD: multiple parameter combinations can give the
exact same MSD (see Beauchemin et al
(2007), Textor
et al (2013)). This means we cannot fit all three parameters based
on the MSD. Here, we’ll set \(t_\text{pause}\) to the literature value of
0.5 minute, and fit the other parameters \(v_\text{free}\) and \(t_\text{free}\) based on the MSD to show
how this works:
# analytical formula
beauchemin.msd <- function( x, M=60, t.free=2, t.pause=0.5, dim=2 ){
v.free <- sqrt( 6 * M * (t.free+t.pause) ) / t.free
M <- (v.free^2 * t.free^2) / 6 / (t.free + t.pause)
multiplier <- rep(1/3, length(x))
xg <- x[x<t.free]
multiplier[x<t.free] <- 1/3 * (xg/t.free)^3 - (xg/t.free)^2 + (xg/t.free)
2 * M * x * dim - 2 * M * dim * t.free * multiplier
}
# Again, fit on the first fitThreshold min, before confinement becomes an issue.
beaucheminFit <- nls( mean ~ beauchemin.msd( dt, M, t.free ),
msdData[ msdData$dt < fitThreshold, ],
start=list(M=30, t.free=1))
# extract parameters M and t.free
M <- coef(beaucheminFit)[["M"]] # in um^2/min
t.free <- coef(beaucheminFit)[["t.free"]] # in min
If everything went correctly, \(M\)
should not be too different from the diffusion coefficient \(D\) we fitted earlier:
## D M
## 19.91724 19.84664
Indeed, they are very similar, which is good news.
Using these fitted values, get the MSD and compare to data:
# Get the fitted MSD using these values
msdBeauchemin <- data.frame(
dt = msdBrownian$dt,
mean = beauchemin.msd( msdBrownian$dt, M, t.free )
)
# plot
ggplot( msdBeauchemin, aes( x = dt, y = mean) ) +
geom_point( data = msdData, shape = 21, aes( fill = dt < fitThreshold ), color = "blue", show.legend = FALSE ) +
geom_line( )+
scale_fill_manual( values = c( "TRUE" = "blue", "FALSE" = "transparent" ) ) +
scale_x_log10( expand = c(0,0) ) +
scale_y_log10( expand = c(0,0) ) +
labs( color = NULL, x = expression( Delta*"t (min)"),
y = expression( "displacement"^2*"("*mu*"m"^2*")") , fill = NULL ) +
theme_bw()
This model seems to fit the data much better!
Now, let’s use the discovered parameters to simulate some tracks on a
longer time scale:
# t.free has been fitted; set t.pause to its fixed value and
# compute v.free from the fitted M:
t.pause <- 0.5 # fixed
v.free <- sqrt( 6 * M * (t.free+t.pause) ) / t.free
# simulate 3 tracks using the fitted parameters.
# we don't use the arguments p.persist, p.bias, bias.dir, and taxis.mode,
# since these are not parameters of the original Beauchemin model.
tau <- timeStep( Txy )
beauchemin.tracks <- simulateTracks( 3,
beaucheminTrack( sim.time = 10*max( timePoints( Txy) ),
delta.t = tau,
t.free = t.free,
v.free = v.free,
t.pause = t.pause ) )
plot( beauchemin.tracks, main = "Beauchemin tracks")
