Single cell motility

Model setup

The model definition file Demo_CellMotility.xml initializes a 2D hexagonal lattice of Lx=200 and Ly=200 with a single cell of target volume 200 \(nodelength^2\) at starting position x=100, y=100 and initial orientation (defined by PropertyVector dir) along the x-axis (dir.phi=0). From t=0...500 with dt=1 per Monte Carlo lattice update of the CPM, this cell performs a persistent random walk imposed by the DirectedMotion plugin with parameters Global/motion_strength and PropertyVector dir. The angle of dir is changed every dt=1 by a normally distributed random variable with mean = 0 and standard deviation given by the parameter Global/noise_level. For analysis purposes, a MotilityReporter calculates finite differences (dt=1) of cell center coordinates along the evolving trajectory and provides a PropertyVector velocity. The Logger Plugin writes a text file with columns time,cell.id,cell.center.x,cell.center.y,velocity.x,velocity.y,velocity.abs and one row only every 100 time steps to emulate the coarser time resolution of image acquisition workflows. On a single CPU core this model costs approximately 0.1sec total wall clock time (plots being disabled). The two parameters Global/motion_strength and Global/noise_level control the statistics of the cell’s persistent random walk.

Ground truth data For the parameter choice Global/motion_strength=1.0 and Global/noise_level=0.1, N=100 repetitions for different random number sequences have been simulated and mean square displacement (MSD) and direction autocorrelation functions (DAC) have been computed from the 100fold coarser time sampling of the simulated cell trajectories. Note, these summary statistics are not (yet) computed within Morpheus but post-hoc in Python and manually entered as np.array measured_data. For plots of these data, see the bottom of the Python notebook.

Parameter recovery study

Using FitMultiCell and adapting the workflow from MinimalExample.ipynb, the parameters are recovered from uniform priors over both parameter intervals from 0.1 to 10 fold the true values. There is some sign of parameter correlations stemming from the effective diffusion coefficient in the MSD statistics such that more random trajectories need a higher speed. To obtain MSD and DAC summary statistics, an own distance measure is computed by the FFT algorithm from the tidynamics module for each single cell trajectory (N=1 and not an ensemble of N=100 as for the ground truth data) with sliding time windows over the 5 temporal data points (tau=0,100,200,300,400 over pairs of data from the trajectory sampled at t=0,100,200,300,400,500 after discarding a burn-in interval of the CPM cell shape initialisation till t=100, hence the data table arguments [1:] in the Python notebook). The meta parameters are set to population_size=20 and 7 epochs.

Details of summary statistics

def distanceMSD(val1, val2):
    d = np.sum(np.abs(tidynamics.msd(np.column_stack([val1['IdSumstat__cell.center.x'][1:],val1['IdSumstat__cell.center.y'][1:]]))\
                      - val2['IdSumstat__MSD']))
    return d
def distanceDAC(val1, val2):
    d = np.sum(np.abs(tidynamics.acf(np.column_stack([val1['IdSumstat__velocity.x'][1:]/val1['IdSumstat__velocity.abs'][1:],val1['IdSumstat__velocity.y'][1:]/val1['IdSumstat__velocity.abs'][1:]]))\
                      - val2['IdSumstat__DAC']))
    return d
distance = pyabc.AggregatedDistance([distanceMSD, distanceDAC])

import pyabc
import fitmulticell as fmc
import scipy as sp
from scipy import stats
import numpy as np
import pandas as pd
import os
import tempfile
%matplotlib inline
import matplotlib.pyplot as plt
from pathlib import Path
import tidynamics                     # to get sliding history stats in N*logN instead of N^2

file_ = str((Path(os.getcwd())) / 'models' / 'Demo_CellMotility.xml')
par_map = {'motion_strength': './Global/Constant[@symbol="motion_strength"]',
           'noise_level': './Global/Constant[@symbol="noise_level"]'}
model = fmc.model.MorpheusModel(
    file_, par_map=par_map,
    executable="/usr/local/bin/morpheus",
    show_stdout=False, show_stderr=True,
    raise_on_error=False)

true_pars = {'motion_strength': 1.0, 'noise_level': 0.1}
limits = {key: (0.1 * val, 10 * val) for key, val in true_pars.items()}

# generate ground truth data for recovery study
# simulated_data = model.sample(true_pars)

# experimentally measured data can be entered manually in this format
# the coarse time resolution of 100 times the simulation step size emulates observability limitations of experiments
measured_data = {'time': np.array([0, 100, 200, 300, 400]),\
                 'MSD': np.array([0, 1640, 5330, 10400, 15800]),\
                 'DAC': np.array([1, 0.49, 0.19, 0.07, 0.07])}

prior = pyabc.Distribution(**{key: pyabc.RV("uniform", lb, ub - lb)
                              for key, (lb, ub) in limits.items()})

# manually defined summary statistics

def distanceMSD(val1, val2):
    d = np.sum(np.abs(tidynamics.msd(np.column_stack([val1['cell.center.x'][1:],val1['cell.center.y'][1:]]))\
                      - val2['MSD']))
    return d

def distanceDAC(val1, val2):
    d = np.sum(np.abs(tidynamics.acf(np.column_stack([val1['velocity.x'][1:]/val1['velocity.abs'][1:],val1['velocity.y'][1:]/val1['velocity.abs'][1:]]))\
                      - val2['DAC']))
    return d

distance = pyabc.AggregatedDistance([distanceMSD, distanceDAC])

abc = pyabc.ABCSMC(model, prior, distance, population_size=20)
db_path = "sqlite:///" + os.path.join(tempfile.gettempdir(), "test1.db")
history = abc.new(db_path, measured_data)

ABC.Sampler INFO: Parallelize sampling on 12 processes.
ABC.History INFO: Start <ABCSMC id=10, start_time=2022-01-20 12:36:30>

abc.run(max_nr_populations=7)

ABC INFO: Calibration sample t = -1.
ABC INFO: t: 0, eps: 2.46221755e+04.
ABC INFO: Accepted: 20 / 58 = 3.4483e-01, ESS: 2.0000e+01.
ABC INFO: t: 1, eps: 2.28795148e+04.
ABC INFO: Accepted: 20 / 56 = 3.5714e-01, ESS: 1.8279e+01.
ABC INFO: t: 2, eps: 2.11752622e+04.
ABC INFO: Accepted: 20 / 72 = 2.7778e-01, ESS: 1.5953e+01.
ABC INFO: t: 3, eps: 1.96163649e+04.
ABC INFO: Accepted: 20 / 131 = 1.5267e-01, ESS: 1.2051e+01.
ABC INFO: t: 4, eps: 1.78708199e+04.
ABC INFO: Accepted: 20 / 360 = 5.5556e-02, ESS: 1.8161e+01.
ABC INFO: t: 5, eps: 1.69895975e+04.
ABC INFO: Accepted: 20 / 485 = 4.1237e-02, ESS: 1.4763e+01.
ABC INFO: t: 6, eps: 1.58884327e+04.
ABC INFO: Accepted: 20 / 795 = 2.5157e-02, ESS: 1.6402e+01.
ABC INFO: Stop: Maximum number of generations.
ABC.History INFO: Done <ABCSMC id=10, duration=0:00:24.166020, end_time=2022-01-20 12:36:55>

<pyabc.storage.history.History at 0x7ff8ad297a10>

h = pyabc.History(db_path)
pyabc.visualization.plot_epsilons(h)
df, w = h.get_distribution(t=h.max_t)
pyabc.visualization.plot_kde_matrix(df, w, limits=limits, refval=true_pars)
plt.show()

# summary statistics for N repetitions of the model simulation and plotting

fit_pars = {'motion_strength': 1.0, 'noise_level': 0.1}
N = 100 # runs to average over, make sure random_seed is off in Morpheus
time = np.arange(5)
meantival1 = np.zeros(5)
meantival2 = np.zeros(5)
semtival1 = np.zeros(5)
semtival2 = np.zeros(5)
count = 0

for i in range(N):
    simulated_data = model.sample(fit_pars)
    val1vec = np.column_stack([simulated_data['cell.center.x'][1:],simulated_data['cell.center.y'][1:]])
    meantival1 += tidynamics.msd(val1vec)
    semtival1 += (tidynamics.msd(val1vec))**2
    val2vec = np.column_stack([simulated_data['velocity.x'][1:]/simulated_data['velocity.abs'][1:],simulated_data['velocity.y'][1:]/simulated_data['velocity.abs'][1:]])
    meantival2 += tidynamics.acf(val2vec)
    semtival2 += (tidynamics.acf(val2vec))**2
    count += 1

meantival1 /= count
semtival1 = np.sqrt(((semtival1 / count) - (meantival1)**2) / count)
meantival2 /= count
semtival2 = np.sqrt(((semtival2 / count) - (meantival2)**2) / count)

print('mean MSD = ',meantival1)
print( 'sem MSD = ',semtival1)
print('mean DAC = ',meantival2)
print( 'sem DAC = ',semtival2)

# plot MSD
plt.plot(measured_data['time'], meantival1, label='Simulation mean')
plt.legend()
plt.title('mean square displacement')
plt.xlabel(r'$\tau$')
plt.ylabel(r'mean square displacement')
plt.show()

# plot DAC
plt.plot(measured_data['time'], meantival2, label='Simulation')
plt.legend()
plt.title('direction autocorrelation')
plt.xlabel(r'$\tau$')
plt.ylabel(r'$\langle v(t) v(t+\tau) \rangle$')
plt.show()

mean MSD =  [-2.56113708e-11  1.55992635e+02  5.99555390e+02  1.34131782e+03
  2.43575981e+03]
sem MSD =  [2.88935230e-12 2.22951285e+00 1.00462531e+01 2.27205059e+01
 3.43561491e+01]
mean DAC =  [1.         0.14711839 0.10313133 0.02550863 0.31417727]
sem DAC =  [6.14390615e-09 3.53107323e-02 5.32866381e-02 6.12744495e-02
 5.66760041e-02]