In this chapter, we will load the pcdl library and use its TimeStep class to load the data snapshot output/00000012, from the 2D time series test dataset (https://github.com/elmbeech/physicelldataloader/blob/master/output_2d.tar.gz).

First, please install the latest version of physicell data loader (pcdl), as described in the HowTo chapter.

And, if not already done so, have a quick read through the pcdl background information.

To runs this tutorial, you can either work with the data that is currently in your output folder, or you can install the 2D unit test dataset into your PhysiCell output folder, by executing the following command sequence.

⚠ Warning: if you run this sequence, all data currently in your PhysiCell/output folder will be overwritten!

cd path/to/PhysiCell

make data-cleanup
python3 -c"import pathlib, pcdl, shutil; pcdl.install_data(); s_ipath=str(pathlib.Path(pcdl.__file__).parent.resolve()/'output_2d'); shutil.copytree(s_ipath, 'output', dirs_exist_ok=True)"

By default, all data related to the snapshot is loaded. For speed and less memory usage, it is however possible to only load the essential (output xml and cell mat data), and exclude microenvironment data, graph data, and possible PhysiBoss data.
For custom_data variables it is possible to specify data types, apart from the generic float type, namely: int, bool, and str. This can be done too, after the data is loaded, using the mcds.custom_data_astype function.
For paths, in general, unix (slash) and windows (backslash) notation will work.

The basic way to load a mcds time step object:

import pcdl  # load the physicell data loader library
print('pcdl version:', pcdl.__version__)  # it is easy to figure out which pcdl version you run

mcds = pcdl.TimeStep('output/output00000012.xml')  # loads the whole snapshot: the xml and all related mat and graph files

The fine tuned way of loading a mcds time step object:

import pcdl  # load the physicell data loader library
print('pcdl version:', pcdl.__version__)  # it is easy to figure out which pcdl version you run

mcds = pcdl.TimeStep(s_pathfile, custom_data_type={}, microenv=True, graph=True, physiboss=True, settingxml='PhysiCell_settings.xml', verbose=True)

The verbosity for pcdl shell output you can tune, even after loading the data.

mcds.set_verbose_false()

mcds.set_verbose_true()

All loaded data lives in mcds.data dictionary.
As in the original python-loader, we tried to keep everything organized inside this dictionary.
Regarding the original python-loader, the structure has slightly changed.\

In pcdl, all data is accessible by functions. There should be no need to fetch data directly from the mcds.data dictionary!

Anyhow, let's take a look at what we actually have in here.

# main data branches
sorted(mcds.data.keys())  # metadata, mesh, substrate (microenvironment), cell (agent), and raw_substrate and raw_cell data.

# metadata
sorted(mcds.data['metadata'].keys())  # multicellds_version, physicell_version, simulation time (current_time), runtime (current_runtime), time stamp (created), time unit ('time_units', 'runtime_units'), spatial unit ('spatial_unit'), and all units specified (ds_unit).

# mesh
sorted(mcds.data['mesh'].keys())  # voxel (ijk), mesh (mnp), and position (xyz) range, axis, spacing, voxel volume, coordinate, and grid.

# microenvironment
sorted(mcds.data['substrate'].keys())  # df_conc, df_substarte, ds_substrate, ls_substarte
print(mcds.data['substrate']['ls_substarte'] # list of all processed substrates
print(mcds.data['substrate']['ds_substrate'] # dictionary with substrate ID label mapping
print(mcds.data['substrate']['df_substarte'] # pandas data frame with substrate decay_rate and diffusion_coefficient.
print(mcds.data['substrate']['df_conc'] # pandas data frame with the actual substrate concentrations in each voxel.

# cell
sorted(mcds.data['cell'].keys())  # dei_graph, df_cell, ds_celltype, ls_cellattr, ls_celltype
print(mcds.data['cell']['ls_celltype'])  # list of all processed cell types
print(mcds.data['cell']['ls_cellattr'])  # list of all outputted (df_cell) cell type related attributes
print(mcds.data['cell']['ds_celltype'])  # dictionary with cell type ID label mapping
print(mcds.data['cell']['df_cell'])  # pandas data frame with each cell's attribute values
sorted(mcds.data['cell']['dei_graph'].keys())  # dictionary with neighbor_cells, attached_cells, and spring_attached_cells graph information.

Once again, loud, for the ones in the back, in pcdl >= version 3, all data is accessible by functions. There should be no need to fetch data directly from the mcds.data dictionaries!

We will explore these functions in the upcoming sections.

Fetch the data's MultiCellDS version, and the PhysiCell version the data was generated.

mcds.get_multicellds_version()  # will return a string like MultiCellDS_2 or MultiCellDS_0.5

mcds.get_pcdl_version()  # will return a string like pcdl_4.1.3

mcds.get_physicell_version()  # will return a string like PhysiCell_1.10.4 or BioFVM_1.1.7

Fetch simulation time, runtime, and time stamp when the data was processed.

mcds.get_time()   # will return a float value like 720.0

mcds.get_runtime()  # will return a float value like 15.596373

mcds.get_timestamp()  # will return a sting like 2022-10-19T01:12:01Z

Finally, it is possible to retrieve a dictionary that lists all units from all tracked variables, from metadata, mesh, substrate, and cells.

mcds.get_unit_dict()  # will return a dictionary, which maps the variables to the units specified in the setting.xml.

# sorted alphabetically
sorted(mcds.get_unit_dict().items())

In the loaded dataset, only one substrate, oxygen, was part of the simulation.
However, let's have a look at the substrate microenvironment related functions.

We can retrieve a list of all substrates processed in the simulation, ordered by substrate ID.

mcds.get_substrate_list()  # ['oxygen']

We can retrieve a dictionary that maps substrate indexes to labels.

mcds.get_substrate_dict()  # {'0': 'oxygen'}

# substrates sorted alpahbetically
sorted(mcds.get_substrate_dict().values())  # ['oxygen']

We can also retrieve a pandas dataframe with the parameters (decay_rate, diffusion_coefficient) that were set for each substrate.

df = mcds.get_substrate_df()
df.head()

Pandas mimics the R computer language, a computer language written for statistical computing. As such, pandas provides us with the spreadsheet like DataFrame and the Series data type.

For analyzing the microenvironment statistically, a pandas dataframe with mesh center coordinates (mnp), voxel coordinates (ijk), and all substrate concentration values can be retrieved.

df_conc = mcds.get_conc_df()
df_conc.info()

Substrates that have in the whole domain the same concentration, carry for that time step no information. Let's filter for substrates that have at least 2 different values over the whole domain.

df_conc = mcds.get_conc_df(values=2)
df_conc.info()  # oxygen concentration varies over the domain

Let's filter for every substrate concentrations in voxel i == 2, j == 1, and k == 0.

df_conc.loc[(df_conc.voxel_i == 2) & (df_conc.voxel_j == 1) & (df_conc.voxel_k == 0), :]

Please have a look at TUTORIAL_python3_pandas.md to learn more.

For substrate concentration visualization, matplotlib contour and contourf plots, for any substrate, through any z_slice can be retrieved.
The mcds.plot_contour function has many parameters to fine tune the plot. Please have a look at its docstring to learn more.

fig = mcds.plot_contour('oxygen')
fig.show()

help(mcds.plot_contour)

The mcds.plot_contour function output can combine with the mcds.plot_scatter output.
Please have a look at TUTORIAL_python3_matplotlib.md to learn more.

For substrate concentration visualization in 3D a rectilinear grid vtk file can be retrieved.
This file can be analyzed, for example with the Paraview software.

mcds.make_conc_vtk()

Please have a look at pcdl and python3 and vtk and TUTORIAL_paraview.md to learn more.

In the loaded dataset, only one cell type, cancer_cell, was part of the simulation.
However, let's have a look at the cell agent related functions.

We can retrieve a list of all cell types processed in the simulation, ordered by cell type ID.

mcds.get_celltype_list()  # ['default']

We can retrieve a dictionary that maps cell type IDs to labels.

mcds.get_celltype_dict()  # {'0': 'default'}

And we can retrieve an alphabetically ordered list with all cell agent attributes outputted by PhysiCell, loadable with mcds.get_cell_df().

mcds.get_cell_attribute_list()  # ['apoptotic_phagocytosis_rate', ...]

Pandas mimics the R computer language, a computer language written for statistical computing. As such, pandas provides us with the spreadsheet like DataFrame and the Series data type.

For analyzing the cell agent data statistically, a pandas dataframe with the agent's position (xyz), the related mesh center coordinate (mnp), the related voxel coordinate (ijk), values from a bunch of cell attributes, the cell surrounding substrate concentrations, and some additional attributes like the voxel's cell density and the xyz vector length, and, if applicable, physiboss output, can be retrieved.

df_cell = mcds.get_cell_df()
df_cell.info()

df_cell.shape  # (992, 95)  this means: 992 agents, 96 tracked variables

Cell attributes, that carry in all agents the same value, carry no information. Let's filter for variables that carry at least 2 different values.

df_cell = mcds.get_cell_df(values=2)
df_cell.shape # (992, 40) this means: 992 agents, 40 tracked variables

Let's filter for cells in voxel i2 j1 k0.

df_cell.loc[(df_cell.voxel_i == 2) & (df_cell.voxel_j == 1) & (df_cell.voxel_k == 0), :]  # cells: 5, 7, 39

Please have a look at TUTORIAL_python3_pandas.md to learn more.

To be able to analyze cell agent data the same way as single cell RNA seq data is analyzed, pcdl has a function to translate cell agent data into anndata format. Anndata is the de facto standard for sc RNA seq data analysis in python. Anndata is the backbone of the scverse (single cell universe) project.

ann = mcds.get_anndata(values=2)
print(ann)  # AnnData object with n_obs × n_vars = 81 × 50
            #     obs: 'z_layer', 'time', 'cell_type', 'chemotaxis_index', 'current_phase', 'cycle_model'
            #     uns: 'neighbor'
            #     obsm: 'spatial'
            #     obsp: 'physicell_neighbor_conectivities', 'physicell_neighbor_distances'

variables:

ann.var_names  # numerical cell attributes:  Index(['attachment_rate', ...], dtype='object')

observation:

ann.obs_names  # cell IDs: Index(['10', ...], dtype='object', name='ID', length=992)

ann.obs_keys()  # categorical cell attributes:  ['z_layer', 'time', 'cell_type', 'chemotaxis_index', 'current_phase', 'cycle_model']

ann.obsm_keys()  # cell coordinates: ['spatial']

ann.obsm['spatial']  # the coordinate values (position\_x, position\_y) of the loaded data.

ann.obsp  # cell neighborhood graph data.

unstructured data:

ann.uns_keys()  # ['neighbor']

ann.uns['neighbor']  # metadata about the neighborhood graph.

The output tells us that we have loaded a time step with 81 cell agents and 50 numerical attributes (vars). Further, we have 6 categorical cell agent attributes (obs). We have each cell agent's xy spatial coordinate information (obsm). And we have cell neighbor graph infromation (obsp, uns).

Please have a look at TUTORIAL_python3_scverse.md to learn more.

For cell agent visualization matplotlib scatter plots, for any cell attribute, through any z_slice can be retrieved.
The mcds.plot_scatter function has many parameters to fine tune the plot. Please have a look at its docstring to learn more.

fig = mcds.plot_scatter()
fig.show()

fig = mcds.plot_scatter('pressure')
fig.show()

help(mcds.plot_scatter)

The mcds.plot_contour function output can combine with the mcds.plot_contour output.
Please have a look at TUTORIAL_python3_matplotlib.md to learn more.

For cell agent visualization in 3D glyph vtk files can be retrieved.
This files can be analyzed, for example with the Paraview software.

mcds.make_cell_vtk()

help(mcds.make_cell_vtk)

Please have a look at pcdl and python3 and vtk and TUTORIAL_paraview.md to learn more.

Fore each cell, PhysiCell tracks the cell neighborhood (if cells touch each other, if cells are attached to each other).

All this information can be represented as a graph (as in graph theory). Pcld has a function to dump this information as gml (graph model language) files. This file format can be read by graph libraries, like network and igraph, for downstream analysis.

Cell neighbor touching graph

mcds.make_graph_gml('neighbor')

Cell neighbor attached graph

mcds.make_graph_gml('attached')

Cell neighbor spring attached graph

mcds.make_graph_gml('spring')

Furthermore, there are functions to recall the graph dictionary for each of the tracked graphs.

mcds.get_neighbor_graph_dict()

mcds.get_attached_graph_dict()

mcds.get_spring_graph_dict()

Please have a look at TUTORIAL_python3_graph.md to learn more.

For bioinformatics, substrate and cell agent data can be saved in the University of Oxford's multi spatial analysis MuSpAn or scverse's SpatiaData format.

do_domain = mcds.get_muspan()

sdata = mcds.get_spatialdata()

Please have a look at TUTORIAL_python3_muspan.md TUTORIAL_python3_scverse.md to learn more.

For substrate and cell agent visualization, data can be saved in open microscopy's ome.tiff file format.

For cell agents, the default cell_attribute outputted is the cell ID + 1, results in segmentation masks, although, any numerical (bool, int, float) cell_attribute can be outputted.

mcds.make_ome_tiff()  # mark by cell ID + 1.

mcds.make_ome_tiff('dead')  # mark dead and alive cells.

mcds.render_neuroglancer(mcds.make_ome_tiff())  # load generated ome.tiff files straigth into the WebGL-based Neuroglancer software.

help(mcds.make_ome_tiff)

help(mcds.render_neuroglancer)

The tiff and ome.tiff files can be loaded back in to python as numpy arrays. Besides that, ome.tiff files enables us to study PhysiCell output the same way as commonly fluorescent microscopy data is analyzed by wetlab scientists. Please have a look at TUTORIAL_python3_ometiff.md, TUTORIAL_python3_napari.md, TUTORIAL_fiji_imagej.md, and [TUTORIAL_neuroglancer.md]((https://github.com/elmbeech/physicelldataloader/blob/master/man/TUTORIAL_neuroglancer.md) to learn more.

For data analysis, the functions related to the mesh are most probably the least one you have to deal with. However, these functions are very useful for tool developers.
Let's have a look at these functions anyway!

It is common, but not necessary, that the voxel's width, height, and depth is the same. In fact, in this test dataset, you will find that this is not the case. In the related PhysiCell_settings.xml file, the voxel was specified as 30[nm] high, 20[nm] wide, and 10[nm] deep.

mcds.get_mesh_spacing()  # the same as get_voxel_spacing.

mcds.get_voxel_spacing() # [30.0, 20.0, 10.0]

mcds.get_voxel_volume()  # 30[nm] * 20[nm] * 10[nm] = 6000[nm\*\*3]

Since version 3, for clarity, coordinate labels are distinct for cell position, mesh center, and voxel index:

• x,y,z: stand for cell position coordinates, and are real values.
• m,n,p: stand for mesh center coordinates, and are real values.
• i,j,k: stand for voxel index coordinates, and are unsigned integer values.

For a better understanding, let's fetch start and stop range values for each coordinate type.
Unlike in the python range function, not only the start value, but also the stop value is inclusive.

mcds.get_voxel_ijk_range()  # [(0, 10), (0, 10), (0, 0)]

mcds.get_mesh_mnp_range()  # [(-15.0, 285.0), (-10.0, 190.0), (0.0, 0.0)]

mcds.get_xyz_range()  # [(-30.0, 300.0), (-20.0, 200.0), (-5.0, 5.0)]

• voxel index values stretch from 0 to 10 longitude (i), 0 to 10 latitude (j), the 2D domain is only 1 voxel deep (k).
• mesh center values stretch from -15[nm] to 285[nm] longitude (m), -10[nm] to 190[nm] latitude (n), and the depth mesh center (p) is at zero.
• cells can hold a positioned between -30[nm] and 300[nm] longitude (x), -20[nm] and 200[nm] latitude (y), and -5[nm] and 5[nm] depth (z).

For voxel and mesh centers, we can fetch the axis tick lists.
Each of this function will return a list of 3 numpy arrays, ordered by ijk or mnp, respective.

mcds.get_voxel_ijk_axis()

mcds.get_mesh_mnp_axis()

There are three helper functions for interacting with the coordinate systems: One to figure out if a particular xyz coordinate is still in side the mesh, the two other ones to translate an xyz coordinate into ijk voxel or mnp mesh center indices, respective. The translation functions will by default checks if the given coordinate is in the mesh.

Is the given xyz coordinate in the mesh?

mcds.is_in_mesh(x=0, y=0, z=0)  # True

mcds.is_in_mesh(x=111, y=22, z=-5)  # True

mcds.is_in_mesh(x=111, y=22, z=-5.1)  # False and Warning @ TimeStep.is_in_mesh : z = -5.1 out of bounds: z-range is (-5.0, 5.0)

Translate a xyz position coordinate into a ijk voxel coordinate:

mcds.get_voxel_ijk(x=0, y=0, z=0)  # [0,0,0]

mcds.get_voxel_ijk(x=111, y=22, z=-5)  # [4, 2, 0]

mcds.get_voxel_ijk(x=111, y=22, z=-5.1)  # None and Warning @ TimeStep.is_in_mesh : z = -5.1 out of bounds: z-range is (-5.0, 5.0)

Translate a xyz position coordinate into a mnp mesh center coordinate:

mcds.get_mesh_mnp(x=0, y=0, z=0)  # [0,0,0]

mcds.get_mesh_mnp(x=111, y=22, z=-5)  # [4, 2, 0]

mcds.get_mesh_mnp(x=111, y=22, z=-5.1)  # None and Warning @ TimeStep.is_in_mesh : z = -5.1 out of bounds: z-range is (-5.0, 5.0)

We can even retrieve all mnp mesh center coordinate triplets, and the meshgrids in 2D or 3D shape itself.

The coordinate triplets are returned as numpy array.

# all mnp coordinates
mnp_coordinats = mcds.get_mesh_coordinate()  # numpy array with shape (3, 121)
mnp_coordinats.shape  # (3, 121)

The meshgrids are returned as a 3 way tensor (2D) or 4 way tensor (3D) numpy array.
As shown below, these tensors can be subset in the common numpy way.

2D mesh grid:

mn_meshgrid = mcds.get_mesh_2D()  # numpy array with shape (2, 11, 11)
mn_meshgrid.shape # (2, 11, 11)

m_meshgrid = mn_meshgrid[0]  # or explicit mn_meshgrid[0,:,:]
m_meshgrid.shape  # (11, 11)

n_meshgrid = mn_meshgrid[1]  # or explicit mn_meshgrid[1,:,:]
n_meshgrid.shape  # (11, 11)

3D mesh grid:

mnp_meshgrid = mcds.get_mesh()  # numpy array with shape (3, 11, 11, 1)
mnp_meshgrid.shape # (3, 11, 11, 1)

m_meshgrid = mnp_meshgrid[0]  # or explicit mnp_meshgrid[0,:,:,:]
m_meshgrid.shape  # (11, 11, 1)

n_meshgrid = mnp_meshgrid[1]  # or explicit mnp_meshgrid[1,:,:,:]
n_meshgrid.shape  # (11, 11, 1)

p_meshgrid = mnp_meshgrid[2]  # or explicit mnp_meshgrid[2,:,:,:]
p_meshgrid.shape  # (11, 11, 1)

After you are done checking out the 2D unit test dataset, you can uninstall the datasets and remove the data in the output folder, by executing the following command sequence.

python3 -c"import pcdl; pcdl.uninstall_data()"
make data-cleanup

That's it!