3D toy model

• generate a segmented 3d image of cells
• modify image to create explicit membrane compartments
• combine cells, combine membranes
• create a simple sme model using this geometry
• do an example simulation

Navigation

• Press Space to show the next page
• Press Shift+Space to show the previous page
• Press Escape to zoom out

Utility functions

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation
from mpl_toolkits.mplot3d import Axes3D
import imageio.v3 as iio
import tifffile
import skimage
import sme
from IPython.display import HTML
from scipy import ndimage as ndi

plt.rcParams["figure.figsize"] = (8, 8)

def plot_voxels(ax, img, filled, title):
    # plot voxels from Z,Y,X,rgb image array to matplotlib Axes3D object ax
    # filled is a Z,Y,X boolean array, if negative the corresponding voxel is not displayed
    # ax.voxels assumes x,y,z order, so first swap the x & z axes
    img_voxels = np.swapaxes(img, 0, 2)
    filled_voxels = np.swapaxes(filled, 0, 2)
    ax.voxels(filled=filled_voxels, facecolors=img_voxels, shade=True)
    ax.set_title(f"{title} [voxels]")
    ax.set_xlabel("x")
    ax.set_ylabel("y")
    ax.set_zlabel("z")

def plot_image_stack(ax, img, title):
    # plot Z,Y,X,rgb image array to matplotlib Axes3D object ax as 2d image for each z-slice
    Z, Y, X, nc = img.shape
    x, y = np.meshgrid(np.arange(X), np.arange(Y))
    for i in range(Z):
        ax.plot_surface(
            x,
            y,
            i * np.ones((Y, X)),
            facecolors=img[i],
            rstride=1,
            cstride=1,
            alpha=1.0,
            shade=True,
            linewidth=0,
        )
    ax.set_title(f"{title} [image stack]")
    ax.set_xlabel("x")
    ax.set_ylabel("y")
    ax.set_zlabel("z")

def show_image(image, title="", ignore_index=None):
    # Expects either Z,Y,X array of voxel indices, or Z,Y,X,color array of (r,g,b) voxel colors
    if image.ndim == 3:
        # convert indexed image -> rgb
        img = skimage.color.label2rgb(image)
        # exclude voxels with index index_to_ignore
        filled = image != ignore_index
    else:
        # image already has rgb values for each voxel: ensure they are normalised to 1
        img = image / np.max(image)
        # include all voxels
        filled = np.full(image.shape[0:3], True)
    fig, (ax_v, ax_i) = plt.subplots(1, 2, subplot_kw=dict(projection="3d"))
    fig.tight_layout()
    plot_voxels(ax_v, img, filled, title)
    plot_image_stack(ax_i, img, title)
    plt.show()

def sphere_mask(grid_shape, center, radius, deformation):
    # generate a boolean mask for a sphere with given center, radius and deformation
    Z, Y, X = grid_shape
    z0, y0, x0 = center
    dz, dy, dx = deformation
    z, y, x = np.ogrid[:Z, :Y, :X]
    return dx * (x - x0) ** 2 + dy * (y - y0) ** 2 + dz * (z - z0) ** 2 <= radius**2

def spheres(n_pixels, n_spheres, max_radius, max_deform):
    # generate a segmented image containing n_spheres randomly distributed, sized and deformed
    voxels = np.zeros((n_pixels, n_pixels, n_pixels), dtype=np.uint16)
    for n_sphere in range(1, n_spheres + 1):
        center = np.random.randint(2, n_pixels - 2, 3)
        nuclear_radius = np.random.randint(1, max_radius / 2)
        cell_radius = np.random.randint(1.5 * nuclear_radius, max_radius)
        deformation = np.random.uniform(1 / max_deform, max_deform, 3)
        voxels[sphere_mask(voxels.shape, center, nuclear_radius, deformation)] = (
            n_sphere
        )
    return voxels

Generate segmented input data

• construct a 40x40x40 3d image with 300 randomly distributed, sized and deformed spheres
• each voxel has an index which identifies which sphere (if any) it belongs to

img_indexed = spheres(n_pixels=40, n_spheres=300, max_radius=14, max_deform=1.5)

show_image((img_indexed), "Segmented image", ignore_index=0)

Generate explicit membranes by dilating each cell

• We want to add explicit membrane compartments around each cell.
• To do this we take a mask of each cell individually, dilate it, and select the pixels that differ from the original mask
• Repeating this over all cells and combining the results gives us a mask of membrane compartment pixels

img_membrane_mask = np.zeros(img_indexed.shape).astype(bool)
kernel = ndi.generate_binary_structure(rank=3, connectivity=1)
kernel_size = (3, 3, 3)
kernel = np.ones(kernel_size, dtype=np.uint8)
for index in range(img_indexed.max()):
    img = (img_indexed == index).astype(np.uint8)
    img_membrane_mask |= ndi.binary_dilation(img) != img

show_image(img_membrane_mask, "Membrane voxels", ignore_index=0)

Define cells as any segmented pixel excluding membrane pixels

• Now we select all pixels that were identified as cells
• Then we exclude pixels that are part of the membrane mask to leave a cell mask

img_cell_mask = img_indexed != 0
img_cell_mask = img_cell_mask & (img_cell_mask != img_membrane_mask)

show_image(img_cell_mask, "Cell voxels", ignore_index=0)

Construct segmented geometry image for sme

• From these masks we can construct a segmented geometry image for sme
• Each colour in this image can then be assigned to a compartment in the model
• We export the result as a 3d tiff to be imported into sme

img = np.zeros(img_cell_mask.shape, dtype=np.uint8)
img[img_cell_mask] = 1
img[img_membrane_mask] = 2
tifffile.imwrite("geom3d.tiff", img)

show_image(img, "Geometry image", ignore_index=0)

Create sme model

• This was done using the GUI, starting from the 2d toy model & importing a 3d geometry image
• As in the 2d case: one species in each compartment, intially only non-zero in outside
• Reactions: outside <-> membrane and membrane <-> cell
• Here we open this model, then import the new geometry image generated above

model = sme.Model("3d-toy-model.xml")

model.import_geometry_from_image("geom3d.tiff")

show_image(model.compartment_image, "Model geometry")

Simulate model

• simulate for 60s, storing the results every 30s
• this might take a few minutes

simulation_results = model.simulate(60, 30)

Simulation results

show_image(
    simulation_results[0].concentration_image / 255.0,
    f"Concentrations at t={simulation_results[0].time_point}",
)

show_image(
    simulation_results[1].concentration_image / 255.0,
    f"Concentrations at t={simulation_results[1].time_point}",
)

show_image(
    simulation_results[2].concentration_image / 255.0,
    f"Concentrations at t={simulation_results[2].time_point}",
)