mechanochat.SEM2#

class mechanochat.SEM2(ne_per_cell: int = 20, re: float = 2., rd_ratio: float = 2.5, adata: AnnData | None = None, cluster_key: str = 'leiden', spatial_key: str = 'spatial', xc: ndarray[Any, dtype[_ScalarType_co]] | None = None, ctype: ndarray[Any, dtype[_ScalarType_co]] | None = None, sim_name: str = 'untitled', seed: int = 0)[source]#

Bases: CellBase

Subcellular Element Method 2D

Create a SEM object

Parameters:

  • ne_per_cell (int, default 20) – Number of elements per cell

  • re (float, default 2.0) – Element radius

  • rd_ratio (float, default 2.5) –

    Cell radius-distance ratio

    rd_ratio>2: cell radius < cell distance/2, tissue with gaps

    rd_ratio=2: cell radius = cell distance/2, no gaps (confluent tissue)

    rd_ratio<2: cell radius > cell distance/2, overcrowded

  • adata (Anndata, optional) –

    Anndata with .obsm[spatial_key] for cell coordinates, .obs[cluster_key] for cell types.

    If not provided, xc and ctype are required

  • xc (np.ndarray, optional) – cell coordinates. Ignored, if adata.obsm[spatial_key] is provided

  • ctype (np.ndarray, optional) – cell types. Ignored, if adata.obs[cluster_key] is provided

  • cluster_key (str, default 'leiden') – Key for cell type in .obs

  • spatial_key (str, default 'spatial') – Key for spatial coordinates in .obsm

  • sim_name (str, default 'untitled') – Simulation name

alpha_csr: csr_matrix#

cell-cell adhesion parameter, shape=(nc,nc). Set by get_param()

beta: ndarray[Any, dtype[_ScalarType_co]]#

cell stiffness parameter, shape=(nc,). Set by get_param()

param_info: dict#

scaling coefficients used by get_param() (adh_a/adh_b/sim_a/sim_b/w/stiff_a/stiff_b, percentile caps)

f_cc: csr_matrix#

cell-cell force matrix (symmetric, abs)

he: ndarray[Any, dtype[_ScalarType_co]]#

element potential energy, shape=(ne,). Set by compute_element_energy()

pairs: ndarray[Any, dtype[_ScalarType_co]]#

interaction element pairs [i,j]

pairs[k] = [i,j]

f_pairs: ndarray[Any, dtype[_ScalarType_co]]#

f_pairs[k] stores force j -> i

pairs[k] = [i,j]

dim: int#

spatial dimension

rc: float#

estimated cell radius in spatial units (median Delaunay distance / rd_ratio)

cell_size: ndarray[Any, dtype[_ScalarType_co]]#

relative cell size, shape=(nc,). Capped at 1.25

scale: float#

normalization scaling factor

deltax: ndarray[Any, dtype[_ScalarType_co]]#

normalization scaling offset. shape=(1,dim)

xc: ndarray[Any, dtype[_ScalarType_co]]#

cell coordinates. shape=(nc,dim)

ne_per_cell: ndarray[Any, dtype[_ScalarType_co]]#

number of elements of each cell, shape=(nc,)

ne: int#

element number

ecid: ndarray[Any, dtype[_ScalarType_co]]#

(ne)

Type:

id of cell to which each element belongs. shape

ceidn: ndarray[Any, dtype[_ScalarType_co]]#

first elements id of each cell (indptr). shape=(nc+1)

Elements of cell i can be retrieve by ceidn[i]:ceidn[i+1].

Number of elements of cell i is ceidn[i+1]-ceidn[i]

xe: ndarray[Any, dtype[_ScalarType_co]]#

element coordinates. shape=(ne,dim)

sim_name: str#

simulation name

t: int#

time step

compute_alphashape(alpha: ndarray[Any, dtype[_ScalarType_co]] | float | None = None, ns: int | None = None, r: float | None = None) → None#

Compute alpha-shape for each cell

Add .alphashape: List[AlphaShape]

compute_contact(k: int = 8, d_th: float | None = None, add_key: str = 'contacts') → None#

Compute cell-cell contacts

add .contact_matrix (csr_matrix)

stored in adata.obsp[add_key]

Parameters:

  • k (int, default 8) – Number of neighbors

  • d_th (Optional[float], default None) – Distance threshold

  • add_key (str, default 'contacts') – Key for storing cell-cell contacts matrix to .obsp

compute_distance(method: str, k: int = 3, return_distances: bool = False) → csr_matrix | None#

Compute cell-cell distances in contact matrix

Add .spatial_distances[method] (csr_matrix)

Parameters:

  • method (str) – valid methods: ‘knn’, ‘delaunay’, ‘contact’

  • k (int, default: 3) – k for knn

  • return_distances (bool, default: False) – Whether to return the distances matrix

Returns:

distance_matrix (csr_matrix) – Cell-cell distances if return_distances is True, None otherwise

get_alpha() → ndarray[Any, dtype[_ScalarType_co]] | None#

Get alpha radius of each cell

Returns:

alpha (Union[np.ndarray, None]) – Alpha radius of each cell if alphashapes have been computed, otherwise None

get_area() → ndarray[Any, dtype[_ScalarType_co]] | None#

Get cell areas

Returns:

area (Union[np.ndarray, None]) – Cell areas if alphashapes have been computed, otherwise None

get_box(cid_list, scaling: bool = True, rotation: Dict | None = None) → ndarray[Any, dtype[_ScalarType_co]]#

Bounding box of a group of cells.

Parameters:

  • cid_list (array-like) – Cell indices.

  • scaling (bool, default: True) – Return the box in spatial units (xe * scale + deltax) instead of simulation units.

  • rotation (dict, optional) – Rotate the elements before measuring the box. Keys: ‘center’ (2,) and ‘theta’ (radians).

Returns:

box (numpy.ndarray) – [xmin, ymin, xmax, ymax], in the format expected by select_cell.

get_elements(i: int, scaling: bool = True) → ndarray[Any, dtype[_ScalarType_co]]#

Get elements of cell i

Parameters:

  • i (int) – Cell index

  • scaling (bool)

Returns:

xe (np.ndarray) – Element coordinates of cell i

get_neighbors(center_i: int)#

Indices of the cells in contact with cell center_i.

Requires compute_contact() to have been run.

Parameters:

center_i (int) – Cell index.

Returns:

contact_i (numpy.ndarray) – Indices of the contacting cells.

select_cell(box) → ndarray[Any, dtype[_ScalarType_co]]#

select cell

Parameters:

box ([xl,yl,xu,yu])

Returns:

cid_list (np.ndarray) – array of selected cells’ index

set_color(cmap_name: str | None = None, n_colors: int | None = None) → None#

Set a color map for sem.ctype (cell types)

Parameters:

  • cmap_name (str) – Name of a seaborn color palette (seaborn.color_palette)

  • n_colors (int, optional) – Number of colors to generate. Defaults to the number of cell types.

update_xc() → None#

Update .xc (cell coordinates)

adata: AnnData | None#

Linked annotation data. None if initialized from direct inputs (xc/ctype).

nc: int#

cell number

ctype: ndarray[Any, dtype[_ScalarType_co]]#

cell type code. shape=(nc)

ctype_list: ndarray[Any, dtype[_ScalarType_co]]#

cell type category. shape=(nct)

color_list: ndarray[Any, dtype[_ScalarType_co]]#

colors for cell type category. shape=(nct,3)

contact_matrix: None | csr_matrix#

cell-cell contact matrix. shape=(nc,nc)

spatial_distances: Dict[str, csr_matrix]#

dictionary to store cell-cell distance matrix

alphashape: List[AlphaShape]#

list of alpha-shape object for each cell. len=nc

alphashape_info: Dict[str, bool | float | int | None]#

information of alpha-shape

alpha_radius: ndarray[Any, dtype[_ScalarType_co]]#

alpha radius of each cell

ns_default: int#

number of augment points for alpha-shape

get_param(df_adh: DataFrame, df_stiffness: DataFrame, embedding_key: str = 'X_pca', adh_range: Tuple[float, float] = (0.25, 4.), stiffness_range: Tuple[float, float] = (0.25, 0.5), w: float = 0.75, return_data: bool = False)[source]#

Compute cell mechanical parameters alpha (adhesion) and beta (stiffness) from gene expression data and store them in self.alpha_csr and self.beta.

Computation pipeline

  1. Neighbor detection — build a KD-tree on cell centroids self.xc and find all cell pairs within a Chebyshev distance of 60 (simulation units). Only neighboring pairs contribute to adhesion and similarity terms.

  2. Adhesion score (adhesion_data) — for each ligand-receptor (LR) pair in df_adh, the per-cell ligand activity l and receptor activity r are computed as the product of min-shifted subunit expression levels (supports multi-subunit ligands/receptors encoded as 'gene1_gene2' in the table). The symmetric interaction score for neighbor pair (i, j) is:

    LR_ij = r_i * l_j + l_i * r_j

    Scores from all LR pairs are summed, capped at the 95th percentile to reduce outlier influence, then scaled to adh_range.

  3. Gene expression similarity (sim_data) — cosine similarity between cell embedding vectors (default: PCA, adata.obsm[embedding_key]) for each neighbor pair, computed as:

    cos_sim_ij = (X_i · X_j) / (||X_i|| * ||X_j||)

    Values are scaled to adh_range.

  4. Alpha matrix — the cell-cell adhesion parameter for pair (i, j) is a weighted combination of expression similarity and LR-based adhesion:

    alpha_ij = sim_data_s * (1 - a) + adhesion_data * a

    Stored as a symmetric sparse matrix self.alpha_csr of shape (nc, nc).

  5. Beta vector — mean expression of cytoskeleton-related genes in df_stiffness per cell, capped at the 95th percentile and scaled to stiffness_range. Stored as self.beta of shape (nc,).

Parameters:

  • df_adh (DataFrame) – Ligand-receptor pair table. Column 0 = ligand gene(s), column 1 = receptor gene(s). Multi-subunit complexes should be joined with '_' (e.g. 'ITGA1_ITGB1').

  • df_stiffness (DataFrame) – Cytoskeleton/stiffness gene table. Column 0 contains gene names whose mean expression determines per-cell stiffness.

  • embedding_key (str, optional) – Key in adata.obsm for the cell embedding used to compute gene expression similarity. Default 'X_pca'.

  • adh_range (tuple of float, optional) – (min, max) range that both adhesion_data and sim_data are scaled to before combining. Default (0.25, 4.0).

  • stiffness_range (tuple of float, optional) – (min, max) range that stiffness values are scaled to. Default (0.25, 0.5).

  • w (float, optional) – Mixing weight for adhesion vs. similarity in the alpha parameter. w=1 → pure LR adhesion; w=0 → pure gene expression similarity. Default 0.75.

  • return_data (bool, optional) – If True, return intermediate arrays for inspection. Default False.

Returns:

  • None – Sets self.alpha_csr (sparse (nc, nc) float32) and self.beta ((nc,) float32) in place.

  • row_indices, col_indices, lr, sim_data, stiffness_data (only when ``return_data=True``) – row_indices, col_indices: COO indices of neighbor pairs (symmetric).

    lr: raw (unscaled) symmetric LR adhesion scores.

    sim_data: raw symmetric cosine similarity values.

    stiffness_data: raw stiffness scores.

sim_gpu(rm_intra: float = 2., rm_inter: float = 3., sigma: float = 0.01, dt: float = 0.02, T: int = 1000, enable_tqdm: bool = True) → None[source]#

Implement the SEM simulation on the GPU.

Parameters:

  • rm_intra (float) – Equilibrium distance of the intracellular element-element potential. Sets the packing of the elements inside a cell.

  • rm_inter (float) – Equilibrium distance of the intercellular element-element potential. Also the default contact threshold used by compute_contact and the alpha-shape radius (_get_e_radius).

  • sigma (float) – Amplitude of the thermal noise.

  • dt (float) – Time step.

  • T (int) – Number of steps.

  • enable_tqdm (bool, default: True) – Show the progress bar.

Returns:

  • None

  • Sets the following fields in `sem` (update `.xe (element coordinates)` and .xc (cell centroids),) – advances .t by T, stores rm_intra / rm_inter in .param, and marks the alpha-shapes as outdated.

Notes

Requires get_param() to have been run, which provides .alpha_csr and .beta.

Intercellular interactions are only evaluated for the cell pairs stored in .alpha_csr and, inside such a pair, only for the elements closer than 30 simulation units along each axis.

NUMBA_CUDA_USE_NVIDIA_BINDING=1 must be set before numba is imported,

compute_force(chunk_size: int = 300000000) → None[source]#

Compute the element-element forces of the current configuration.

All element pairs closer than 30 simulation distance (Chebyshev) are enumerated, and the force of each pair is evaluated on the GPU with the same parameters as sim_gpu. Pairs are processed in chunks to bound the GPU memory.

Parameters:

chunk_size (int, default: 300000000) – Number of element pairs sent to the GPU per launch. Lower it if the GPU runs out of memory.

Returns:

  • None

  • Sets the following fields in `sem` –

    • .pairs : (n_pairs, 2) element pair indices, pairs[k] = [i, j]

    • .f_pairs : (n_pairs, dim) the force from j exerted on i for pair k

    • .f_cc and adata.obsp[‘f_cc’]: (n_cell, n_cell) cell-cell mechanical interaction matrix.

Notes

Requires .param.

sim_gpu() must have been run first.

compute_element_energy()[source]#

Compute the Lennard-Jones potential energy of each element.

Sums the intracellular contributions (weighted by beta) and the intercellular ones (weighted by alpha_csr) with the same cutoffs as sim_gpu.

Returns:

  • None

  • Sets the following fields in `sem` –

    • .he (ne,) array of Lennard-Jones potential energy of each element.

compute_cell_forces() → None[source]#

Aggregate the element-element forces into a cell-cell mechanical interaction matrix.

The magnitudes of the forces of all element pairs that belong to two different cells are summed, which gives a symmetric measure of how strongly two cells push against or pull on each other. Entries below 0.1 are dropped.

Returns:

  • None

  • Sets the following fields in `sem` – Adds .f_cc, a symmetric (nc, nc) sparse matrix, and copies it to adata.obsp[‘f_cc’] when the SEM object is linked to an AnnData. Called automatically at the end of compute_force().

f_pairs_tocsr(e_mask: ndarray[Any, dtype[bool_]] | None = None) → Tuple[csr_array, csr_array][source]#

Convert force pairs to CSR sparse matrices. Run .compute_force() first

fx_csr[i, j] = fx_{j->i}

row index i : under force column index j : exert force

Parameters:

e_mask (Optional[NDArray[np.bool_]]) – Boolean mask of shape (ne,) to select subset of elements. If provided, only interactions involving selected elements are included, and the returned matrices have shape (n_selected, n_selected) with re-indexed rows/columns. If None, return full matrices of shape (ne, ne).

Returns:

  • fx_csr (csr_array) – Sparse matrix of x-component forces

  • fy_csr (csr_array) – Sparse matrix of y-component forces

compute_tensor() → None[source]#

compute element stress tensor

run compute_force first

add .e_sigma

compute_pressure() → None[source]#

compute element pressure

run compute_tensor first

add .p

save_sim() → None[source]#

Save the simulation to {sim_name}_{t}.pkl in the working directory.

Stores xe, ceidn, ecid, param, scale and deltax, plus pairs / f_pairs, beta and alpha_csr when they exist. The AnnData object is not saved; restore the object with load_sim() on a SEM built from the same AnnData.

If the file already exists, a warning is issued and _temp is appended to the name instead of overwriting it.

load_sim(sim_name: str, t: int, path: str = '.', rename: bool = True) → None[source]#

Restore a simulation from {path}/{sim_name}_{t}.pkl

Parameters:

  • sim_name (str) – Name of simulation

  • t (int) – Time point

  • path (str, default: '.') – Path to simulation data

  • rename (bool, default: True) – If True, rename the sem to sim_name