ForSys: non-invasive stress inference from time-lapse microscopy

Augusto Borges 1, 2, Jerónimo R. Miranda-Rodríguez 1, 3, Alberto Sebastián Ceccarelli 4, Guilherme Ventura5, Jakub Sedzinski 5, Hernán López-Schier 1,2,6 & Osvaldo Chara 7, 8, 9

1. Unit Sensory Biology, Helmholtz Zentrum München, Munich, Germany

2. Graduate School of Quantitative Biosciences (QBM), Munich, Germany

3. Instituto de Neurobiología, Universidad Nacional Autónoma de México (UNAM), Boulevard Juriquilla 3001, Juriquilla, México

4. Systems Biology Group (SysBio), Institute of Physics of Liquids and Biological Systems (IFLySIB), National Scientific and Technical Research Council (CONICET), University of La Plata, La Plata, Argentina

5. The Novo Nordisk Foundation Center for Stem Cell Medicine (reNEW), University of Copenhagen, Blegdamsvej 3B, 2200, Copenhagen, Denmark

6. Division of Science, New York University Abu Dhabi, Saadiyat Island, United Arab Emirates

7. School of Biosciences, University of Nottingham, Sutton Bonington Campus, Nottingham, LE12 5RD, UK

8. Instituto de Tecnología, Universidad Argentina de la Empresa, Buenos Aires, Argentina

9. Corresponding author: osvaldo.chara@nottingham.ac.uk

Generate an inference in a set of SurfaceEvolver steps

import sys
sys.path.append('..')
import forsys as fs
import os
import matplotlib.pyplot as plt
import numpy as np
import matplotlib as mpl

DATA_FOLDER = os.path.join("data", "in_silico")
RESULTS_FOLDER = os.path.join("results")

max_time = 25

# This is only necessary if you wish to create outputs
# if not os.path.exists(RESULTS_FOLDER):
#     os.makedirs(RESULTS_FOLDER)

Create all the frames by iterating in time.

The real_time is used to get the \(\Delta t\) in the velocity calculation

frames = {}
for time in range(max_time):
    if time == 0:
        real_time = 0
    elif time == 1:
        real_time += 1 * 25000 * 0.005
    else:
        real_time += 1 * 50 * 0.005

    lattice = fs.surface_evolver.SurfaceEvolver(os.path.join(DATA_FOLDER,
                                                             f"step_{time}.dmp"))
    frames[time] = fs.frames.Frame(time,
                                   lattice.vertices,
                                   lattice.edges,
                                   lattice.cells,
                                   time=real_time,
                                   gt=True)

Now create the ForSys object from all the frames. The mesh is created by automatically tracking the vertices through time

forsys = fs.ForSys(frames)

Build and solve the system of equations for the force and pressure

You could solve only the time of interest or all the times with a loop. The b_matrix parameter turns the Dynamic instance on

all_velocities = forsys.get_system_velocity_per_frame()
ave_velocity_to_normalize = np.mean(all_velocities)
for time in range(max_time):
    forsys.build_force_matrix(when=time)
    forsys.solve_stress(when=time,
                        b_matrix="velocity",
                        velocity_normalization=0.07 / ave_velocity_to_normalize,
                        adimensional_velocity=False,
                        method="nnls",
                        use_std=False,
                        allow_negatives=False)
    forsys.build_pressure_matrix(when=time)
    forsys.solve_pressure(when=time, method="lagrange_pressure")

d:\manuscripts\unpublished\protocol_forsys\forsys\.venv\Lib\site-packages\forsys\fmatrix.py:317: UserWarning: Numerically solving due to the following error: Singular matrix
  warnings.warn(f"Numerically solving due to the following error: {e}")

Create system’s plots

fig, axs = plt.subplots(5, 5, figsize=(15, 15))

for time, current_ax in enumerate(axs.reshape(-1)):
    _, ax = fs.plot.plot_inference(forsys.frames[time],
                                    normalized="max",
                                    mirror_y=False,
                                    colorbar=False,
                                    pressure=False,
                                    ax=current_ax)
    current_ax = ax
    current_ax.set_title(f"Time: {time}")
    current_ax.axis("off")


plt.tight_layout()
plt.show()

The pressure can be plotted using the pressure argument

fig, axs = plt.subplots(5, 5, figsize=(15, 15))

for time, current_ax in enumerate(axs.reshape(-1)):
    _, ax = fs.plot.plot_inference(forsys.frames[time],
                                    normalized="max",
                                    mirror_y=False,
                                    colorbar=False,
                                    pressure=True,
                                    ax=current_ax)
    current_ax = ax
    current_ax.set_title(f"Time: {time}")
    current_ax.axis("off")


plt.tight_layout()
plt.show()

Data can be exported through the Frames.get_tensions() method, and compared to the ground truth if it exists

fig, ax = plt.subplots()
time_colormap = plt.get_cmap("cividis")
for time in range(max_time)[::5]:
    stresses = forsys.frames[time].get_tensions()

    plt.scatter(stresses["gt"] / stresses["gt"].mean() , stresses["stress"], color=time_colormap(time / max_time), alpha=0.5)


sm = plt.cm.ScalarMappable(cmap="cividis", norm=plt.Normalize(vmin=0, vmax=max_time))

plt.plot([0, 2], [0, 2], color="black", ls="--")
plt.xlabel("Ground truth")
plt.ylabel("Inferred stress")
plt.colorbar(sm, ax=ax, label="Time")
plt.show()