Extensions to Tyssue: solver updates, topology transitions, and cell behaviour functions

.gitignore

@@ -1,6 +1,6 @@
 *.py[cod]
 
-# C extensions
+# C extensions 
 *.so
 *.lo
 *.la

src/tyssue/behaviors/sheet/basic_events.py

@@ -7,7 +7,8 @@
 import logging
 
 from ...geometry.sheet_geometry import SheetGeometry
-from ...topology.sheet_topology import cell_division
+from ...topology.sheet_topology import cell_division, type1_transition as T1_transition
+from ...topology.base_topology import drop_face
 from ...utils.decorators import face_lookup
 from .actions import (
     decrease,
@@ -248,3 +249,41 @@ def contraction_line_tension(sheet, manager, **kwargs):
         isotropic=True,
         limit=100,
     )
+
+
+def T1Swap(sheet,manager, geom, t1_threshold, multiplier):
+    """
+    A behaviour function of the T1 transition that performs a T1 swap on an edge that is shorter than T1_threshold.
+    """
+    # First, we need to get the joint index over free and east edges,
+    # that is, the indices of edges that is spanning the entire graph without double edges
+    sheet.get_extra_indices()
+    edge_dataframe = sheet.edge_df.loc[sheet.sgle_edges]
+    short_edges = edge_dataframe.loc[(edge_dataframe['length'] < t1_threshold)]
+    # take advantage of the unique ID to help us tracking the edges,
+    # we need to reindex the df in type 1 transition function, otherwise there will be inconsistency between
+    # dataframes and will cause out of index error.
+    unique_list = short_edges['unique_id'].tolist()
+    for ID in unique_list:
+        idx = sheet.idx_lookup(ID,'edge')
+        print(f'Performed T1 swap on edge {idx}')
+        T1_transition(sheet, idx, do_reindex=True, remove_tri_faces=False, multiplier=multiplier)
+        sheet.reset_index()     # removes disconnected vertices and faces
+        geom.update_all(sheet)
+    manager.append(T1Swap, geom= geom, t1_threshold = t1_threshold, multiplier = multiplier)
+
+
+def T2Swap(sheet, manager, crit_area):
+    """
+    A behaviour function of the T2 transition that should be added to the manager during simulation.
+    It removes the face if it is triangular and its area is smaller than crit_area.
+    """
+    face_dataframe = sheet.face_df
+    Face_list = face_dataframe.loc[(face_dataframe['num_sides'] < 4) & (face_dataframe['area'] < crit_area)].index.tolist()
+    for face in Face_list:
+        drop_face(sheet, face)
+    manager.append(T2Swap, crit_area = crit_area)
+
+
+
+

src/tyssue/behaviors/sheet/bilayer_dummy_set.py

@@ -0,0 +1,33 @@
+""" The function in this script auto-controls the edges are active or not for a bilayer tissue sheet."""
+
+def bilayer_dummy_set(sheet):
+    """
+    Set edges as active or dummy for bilayer tissue sheet.
+
+    Parameters
+    ----------
+    sheet : tyssue.Sheet
+        The tissue sheet instance containing face_df and edge_df DataFrames.
+    """
+
+    # Iterate through each edge in the edge DataFrame
+    for i in sheet.edge_df.index:
+        # Check if the edge has an opposite edge (i.e., it's internal)
+        if sheet.edge_df.loc[i, 'opposite'] != -1:
+            # Get the associated cell (face) for this edge
+            associated_cell = sheet.edge_df.loc[i, 'face']
+            # Get the opposite edge index
+            opposite_edge = sheet.edge_df.loc[i, 'opposite']
+            # Get the opposite cell (face) for the opposite edge
+            opposite_cell = sheet.edge_df.loc[opposite_edge, 'face']
+            # If both associated and opposite cells are of class 'STB', set edge as dummy (inactive)
+            if (sheet.face_df.loc[associated_cell, 'cell_class'] == 'STB' and
+                    sheet.face_df.loc[opposite_cell, 'cell_class'] == 'STB'):
+                sheet.edge_df.loc[i, 'is_active'] = 0
+                sheet.edge_df.loc[opposite_edge, 'is_active'] = 0
+            else:
+                # Otherwise, set edge as active
+                sheet.edge_df.loc[i, 'is_active'] = 1
+        else:
+            # For boundary edges, set as active
+            sheet.edge_df.loc[i, 'is_active'] = 1

src/tyssue/behaviors/sheet/cell_activity_events.py

@@ -0,0 +1,125 @@
+"""
+This file contains all behaviour functions that models different cellular activities.
+
+"""
+
+import numpy as np
+import pandas as pd
+from ...geometry.planar_geometry import PlanarGeometry
+from ...topology.sheet_topology import cell_division, type1_transition, remove_face
+
+"""cell proliferation behaviour function
+
+A cell proliferation behaviour function does two thing on the selected cell.
+First, it compares the current cell area with an area threshold to see if the cell large enough for division.
+Secondly: if the cell is large enough, a cell division is performed on the cell; alternatively, if the cell area is 
+smaller than the threshold value, then the target area is added by "growth_rate * dt" to expand the cell.
+"""
+# Note: still needs to improve the function, so it controls cell class change.
+def proliferation(sheet, manager, geom, unique_id, crit_area, growth_rate, dt):
+    idx = sheet.idx_lookup(unique_id,'face') # get the current face index from unique_id.
+    if sheet.face_df.loc[idx, "area"] > crit_area:
+        # restore prefered_area
+        sheet.face_df.loc[idx, "prefered_area"] = 1.0
+        # Do division
+        daughter = cell_division(sheet, idx, geom)
+        # Update the topology
+        sheet.reset_index(order=True)
+        # update geometry
+        geom.update_all(sheet)
+        print(f"cell n°{daughter} is born")
+    else:
+        sheet.face_df.loc[idx, "prefered_area"] *= (1 + dt * growth_rate)
+        manager.append(proliferation, geom = geom, unique_id = unique_id, crit_area = crit_area, growth_rate = growth_rate, dt = dt)
+
+
+"""cell fusion behaviour function
+
+A cell fusion behaviour function is used when a CT is fusing into the STB layer. The cell class of the selection cell 
+should become "STB" at the end of the function.
+First, the edges between the selected cell and its STB neighbours are disabled for edge tension term (coefficient = 0).
+
+Secondly, we need to perform a T1 swap on the edge that connects a boundary vertex and the mutual vertex shared between all 
+STB neighbours and the fusing cell; in this way, the newly fused cell would have an edge that is "open" to "outside".
+
+Thirdly, new dynamic parameters need to be updated to ensure consistent physics rule.
+"""
+def fusion(sheet, manager, geom, unique_id):
+    sheet.get_extra_indices()
+    idx = sheet.idx_lookup('face', unique_id)
+    # store the face index of STB neighbours
+    STB_neighbours = sheet.get_neighbors(idx)
+    STB_neighbours = list(STB_neighbours.intersection(set(sheet.face_df.loc[sheet.face_df['cell_class'] == 'STB'].index)))
+    # Find the edges associated with the fusing face index, filter out the boundary edges.
+    internal_edges = sheet.edge_df[(sheet.edge_df['face'] == idx) & (sheet.edge_df['opposite'] != -1)]
+    # We have to update both the internal arrowed edge and its opposite arrowed edge.
+    for ie in internal_edges.index:
+        sheet.edge_df.loc[ie,'is_active'] = 0
+        sheet.edge_df.loc[sheet.edge_df.loc[ie,'opposite'],'is_active'] = 0
+
+    # Find all the boundary vertices in STB neighbours, based on the opposite == -1 value.
+    STB_boundary_edges = sheet.edge_df[
+        (sheet.edge_df['face'].isin(STB_neighbours)) &
+        (sheet.edge_df['opposite'] == -1)
+        ]
+    STB_boundary_verts = pd.unique(STB_boundary_edges[['srce','trgt']].values.ravel())
+    # Next, extract all the vertices belong to the fusing face, the edge connects a boundary vertex and the fusing face
+    # is the edge that we should perform T1 swap on. We can utilise the variable internal_edges.
+    fusing_face_verts = internal_edges['srce']
+    # We only need to looping over the STB neighbours edges, then find edges connects a boundary vertex to a fusing face vertex.
+    STB_edges = sheet.edge_df[sheet.edge_df['face'].isin(STB_neighbours)]
+    matching_edge = STB_edges[
+        (
+                (STB_edges['srce'].isin(fusing_face_verts)) & (STB_edges['trgt'].isin(STB_boundary_verts))
+        ) |
+        (
+                (STB_edges['trgt'].isin(fusing_face_verts)) & (STB_edges['srce'].isin(STB_boundary_verts))
+        )
+    ]
+    if matching_edge.empty:
+        raise ValueError("No matching edge found between fusing face vertices and STB boundary vertices.")
+    first_edge_index = matching_edge.index[0]
+    New_boundary_edge = type1_transition(sheet,first_edge_index,do_reindex=False, remove_tri_faces=False, multiplier=1.5)
+    # Then make the new boundary edge to be active in tension (was dummy before T1).
+    sheet.edge_df.loc[New_boundary_edge,'is_active'] = 1
+    geom.update_all(sheet)
+
+
+"""cell extrusion behaviour function
+
+A cell extrusion behaviour function has two components.
+
+The first component is to let the STB to detach from the CT layer. This is done by a series of T1 transition on shared
+edges between the selected STB and CTs: Detach STB.
+
+The second component is to let the STB to shed from the layer. The shedding is modelled by cell removal but keep the
+vertices shared between STB units that are still in the system: STB removal.
+"""
+def extrude(sheet, manager, geom, unique_id):
+    idx = sheet.idx_lookup('face', unique_id)
+    remove.face(sheet,idx)
+    geom.update_all(sheet)
+
+def detach(sheet, manager, geom, unique_id):
+    sheet.get_extra_indices()
+    idx = sheet.idx_lookup('face', unique_id)
+    # Identify CT neighbours (non-STB)
+    CT_neighbours = sheet.get_neighbors(idx)
+    CT_neighbours = list(CT_neighbours.intersection(
+        set(sheet.face_df.loc[sheet.face_df['cell_class'] != 'STB'].index)
+    ))
+    # Find internal edges of the detaching face
+    internal_edges = sheet.edge_df[(sheet.edge_df['face'] == idx) & (sheet.edge_df['opposite'] != -1)]
+    # Filter internal edges that are shared with CT neighbours
+    shared_edges = internal_edges[
+        sheet.edge_df.loc[internal_edges['opposite'], 'face'].isin(CT_neighbours).values
+    ]
+    if shared_edges.empty:
+        raise ValueError("No mutual edge found between detaching face and CT neighbours.")
+    # Perform T1 transitions on each shared edge and update geometry
+    for edge_idx in shared_edges.index:
+        new_edge_idx = type1_transition(sheet, edge_idx, do_reindex=False, remove_tri_faces=False, multiplier=1.5)
+        sheet.edge_df.loc[new_edge_idx, 'is_active'] = 1
+        geom.update_all(sheet)
+    # Optionally extrude the detached face
+    manager.append(extrude, geom = geom, unique_id = unique_id)

src/tyssue/behaviors/sheet/cell_class_events.py

@@ -0,0 +1,67 @@
+"""
+Event module for cell class transition rules, modify the details accordingly to your model.
+=======================
+
+"""
+
+import numpy as np
+from ...geometry.planar_geometry import PlanarGeometry
+from ...topology.sheet_topology import cell_division
+
+def cell_cycle_transition(sheet, manager, dt, p_recruit=0.1, G2_duration=0.4, G1_duration=0.11):
+    """
+    Controls cell class state transitions for cell cycle based on timers and probabilities.
+
+    Parameters
+    ----------
+    sheet: tyssue.Sheet
+        The tissue sheet.
+    manager: EventManager
+        The event manager scheduling the behaviour.
+    face_id: Integer
+        ID of the cell being controlled.
+    p_recruit: float
+        Probability for an 'S' cell to be recruited to 'G2'.
+    dt: float
+        Time step increment.
+    G2_duration: float
+        Fixed duration cells stay in G2 phase.
+    G1_duration: float
+        Fixed duration cells stay in G1 phase.
+    """
+    # Generate two df that are cells that need to change its cell class or stay in its current class.
+    cells_stay_in_class = sheet.face_df.loc[sheet.face_df['timer'] > 0].index.tolist()
+    cells_to_change = sheet.face_df.loc[sheet.face_df['timer'] <= 0]
+
+    # For cells that still needs to elapse the timer, just reduce the timer by dt
+    for cell in cells_stay_in_class:
+        sheet.face_df.loc[cell,'timer'] -= dt
+
+    # Then we change the cell type based on the cell cycle diagram.
+    G1_cells = cells_to_change.loc[cells_to_change['cell_class'] == 'G1'].index.tolist()
+    S_cells = cells_to_change.loc[cells_to_change['cell_class'] == 'S'].index.tolist()
+    G2_cells = cells_to_change.loc[cells_to_change['cell_class'] == 'G2'].index.tolist()
+    M_cells = cells_to_change.loc[cells_to_change['cell_class'] == 'M'].index.tolist()
+
+    # For cells in G1_cells indices, we simply change the cell class to S
+    sheet.face_df.loc[G1_cells, 'cell_class'] = 'S'
+
+    # For cells in S_cells , we need to loop and decide based on the random number generated for if it moves into G2.
+    for cell in S_cells:
+        if np.random.rand() < p_recruit:
+            sheet.face_df.loc[cell, 'cell_class'] = 'G2'
+            sheet.face_df.loc[cell, 'timer'] = G2_duration
+
+    # For cells in G2_cells, we move them into M class
+    sheet.face_df.loc[G2_cells, 'cell_class'] = 'M'
+
+    # For cells in M_cells, they need to be undergone cell division
+    for cell in M_cells:
+        daughter = cell_division(sheet, mother=cell, geom = PlanarGeometry )
+        # Set parent and daughter to G1 with G1 timer, note that variable daughter is the index of the new row already.
+        sheet.face_df.loc[cell, 'cell_class'] = 'G1'
+        sheet.face_df.loc[daughter, 'cell_class'] = 'G1'
+        sheet.face_df.loc[cell, 'timer'] = G1_duration
+        sheet.face_df.loc[daughter, 'timer'] = G1_duration
+
+    manager.append(cell_cycle_transition, dt = dt, p_recruit = p_recruit,G2_duration = G2_duration, G1_duration = G1_duration)

src/tyssue/core/objects.py

@@ -570,7 +570,8 @@ def get_orbits(self, center, periph):
         return orbits
 
     def idx_lookup(self, elem_id, element):
-        """returns the current index of the element with the `"id"` column equal to `elem_id`
+        """returns the current index of the element with the `"id"` column equal to `elem_id`;
+        thse stable IDs are called "ID" in face data frame, but called "unique_id" in vertex and edge data frames
 
         Parameters
         ----------
@@ -579,7 +580,10 @@ def idx_lookup(self, elem_id, element):
         element : {"vert"|"edge"|"face"|"cell"}
           the corresponding dataset.
         """
-        df = self.datasets[element]["id"]
+        if element in ['vert', 'edge','face']:
+            df = self.datasets[element]['unique_id']
+        else:
+            df = self.datasets[element]["id"]
         idx = df[df == elem_id].index
         if len(idx):
             return idx[0]

src/tyssue/draw/bilayer_drawing_tool.py

@@ -0,0 +1,40 @@
+# This script contains a function that auto updates the drawing specifications for a bilayer tissue sheet.
+
+import numpy as np
+
+def bilayer_draw_spec_update(sheet, specs):
+    """
+    Update drawing specifications for bilayer tissue sheet.
+
+    Parameters
+    ----------
+    sheet : tyssue.Sheet
+        The tissue sheet instance containing face_df and edge_df DataFrames.
+    specs : dict
+        The current drawing specifications dictionary to be updated.
+    """
+
+    # --- FACE COLOR UPDATE ---
+    # Use NumPy vectorization to assign colors:
+    # If 'cell_class' == 'STB', then color = 0.7
+    # Else, then color = 0.1
+    sheet.face_df['color'] = np.where(
+        sheet.face_df['cell_class'] == 'STB', 0.7, 0.1
+    )
+
+    # Update the specs dictionary with the new face colors
+    specs['face']['color'] = sheet.face_df['color'].to_numpy()
+
+    # Set transparency (alpha) for faces
+    specs['face']['alpha'] = 0.2
+
+    # --- EDGE WIDTH UPDATE ---
+    # Use NumPy vectorization to assign edge widths:
+    # If 'is_active' == 0, then width = 2
+    # Else, then width = 0.5
+    sheet.edge_df['width'] = np.where(
+        sheet.edge_df['is_active'] == 0, 2, 0.5
+    )
+
+    # Update the specs dictionary with the new edge widths
+    specs['edge']['width'] = sheet.edge_df['width'].to_numpy()

src/tyssue/draw/plt_draw.py

@@ -129,7 +129,39 @@ def create_gif(
         plt.close(fig)
 
     try:
-        subprocess.run(["convert", (graph_dir / "movie_*.png").as_posix(), output])
+        # Expand pixels manually (cross-platform safe)
+        pngs = sorted(graph_dir.glob("movie_*.png"))
+        png_paths = [str(p) for p in pngs]
+
+        # Detect the correct ImageMagick executable:
+        # IM6: convert
+        # IM7: magick convert
+        def find_imagemagick():
+            # Try IM7 first
+            try:
+                subprocess.run(["magick", "-version"], stdout=subprocess.DEVNULL, stderr=subprocess.DEVNULL)
+                return ["magick", "convert"]
+            except Exception:
+                pass
+
+            # Try IM6
+            try:
+                subprocess.run(["convert", "-version"], stdout=subprocess.DEVNULL, stderr=subprocess.DEVNULL)
+                return ["convert"]
+            except Exception:
+                pass
+
+            raise RuntimeError("ImageMagick executable not found on this system.")
+
+        im_cmd = find_imagemagick()
+
+        # Run ImageMagick safely
+        try:
+            subprocess.run(im_cmd + png_paths + [output], check=True)
+        except Exception as e:
+            print("Converting didn't work. Make sure ImageMagick is correctly installed.")
+            raise e
+
     except Exception as e:
         print(
             "Converting didn't work, make sure imagemagick is available on your system"