1. Scientific Overview
Three-dimensional multicellular aggregates, commonly classified as organoids, represent an exceptional scientific development in bioengineering and drug screening. Generated in vitro from induced pluripotent stem cells ($iPSCs$) or embryonic stem cells ($ESCs$), these constructs undergo guided self-organization to recapitulate structural features of human organs (e.g., cortical structures, intestinal crypts, nephron networks). Unlike flat two-dimensional monolayers, 3D systems are constrained by diffusion limitations and localized physical mechanical forces that establish complex physiological and biological patterns.
In physical biological matrices, cellular development and lineage commitment are determined by chemical concentration profiles. To model these phenomena, we implement a coupled computational engine combining soft-body physical forces and biochemical laws. Consider a spatial domain where the nutrient concentration (such as dissolved oxygen) $C(\mathbf{x}, t)$ is defined by a standard partial differential equation for reaction-diffusion dynamics:
$$ \frac{\partial C}{\partial t} = D \nabla^2 C - k C $$
where $D$ represents the spatial diffusion coefficient of nutrients within the local cellular extracellular matrix (ECM), and $k$ denotes the metabolic consumption rate of the cell cluster. If an organoid grows beyond the critical diffusion boundary limit ($>200\ \mu\text{m}$), the internal nutrient core concentration falls below a critical threshold $C_{\text{crit}}$, yielding cell cycle arrest and triggering a transition to an apoptotic or necrotic state.
Concurrently, cell division (mitosis) is moderated by contact-mediated mechanics. We model contact inhibition of proliferation, a biological property where physical compression halts replication. As cellular packing density $P = N / \text{Area}$ approaches a saturation parameter $P_{\text{max}}$, the proliferative capacity of pluripotent cells drops:
$$ \Gamma = \Gamma_0 \max\left(0, 1 - \frac{P}{P_{\text{max}}}\right) $$
where $\Gamma_0$ is the unconstrained stem-cell replication velocity. This simulation represents a dynamic interactive agent-based system designed to illustrate these exact biophysical mechanisms.
2. How to Use
The simulator is highly reactive, responding to user-directed spatial physical adjustments and biophysical coefficients. Users can interact through direct visual mechanics and numerical controls:
- Manual Seeding & Excision: Click and drag anywhere inside the black visualizer canvas to seed new stem cell coordinates (green nodes) within the virtual matrix. Right-click and drag (or long-press and drag on touch screens) to trigger spatial physical excision, removing cells from that area.
- Mitogen Injection (+ Mitogen): Triggers cell division among active green stem cells. Mitogen exposure induces mitosis, causing stem cells to divide into neighboring vacant space, provided they are not restricted by local contact inhibition.
- Inducing Differentiation: Forces stochastic lineage commitment, converting pluripotency-retaining stem cells (green) into non-proliferative, stable differentiated somatic cell types (blue).
- Toxin Injection: Introduces chemical-stress agents across the canvas. Randomly targets cells, turning them necrotic (red), which initiates a slow physical structural breakdown and eventual clearance.
- Biophysical Coefficients Sliders: Modify the sliding physical inputs inside the panel to adjust global parameters:
- Fluid Viscosity ($\mu$): Alters mechanical damping. High viscosity yields high friction, slowing cellular displacement.
- Cadherin Cohesion ($k_a$): Sets cell-cell cadherin adhesion strength. High values pull cells into dense, compact spheres, simulating tight cell junctions.
- Core Elasticity ($k_r$): Governs repulsion force when cells overlap, mimicking the elastic cell wall structural response.
- Symmetric Fate Probability ($p_d$): Dictates the stochastic baseline rate of natural, uninduced stem-to-somatic cell differentiation.
- Demo Mode & Baseline Reset: Press "Start Demo" to initiate a dynamic guided tour of physiological scenarios (proliferation, maturation, homeostasis, toxic events, clearance, and recovery). Click "Reset Baseline" to instantly clear the platform and restore default biophysical parameters.
3. Technical Details
Our computational engine resolves agent-based dynamics via modified Verlet integration algorithms coupled to a spatial partitioning grid. Cells are treated as elastomeric soft spheres possessing a finite physical interaction boundary. The position vector of cell agent $i$ is calculated using an implicit velocity-based integration framework:
$$ \mathbf{v}_i(t + \Delta t) = \left( \mathbf{v}_i(t) + \frac{\mathbf{F}_i(t)}{m_i} \Delta t \right) \cdot \mu $$
$$ \mathbf{x}_i(t + \Delta t) = \mathbf{x}_i(t) + \mathbf{v}_i(t + \Delta t) \Delta t $$
where $\mathbf{F}_i(t)$ represents the net physical force vector acting on the cell, $m_i$ is the cellular mass, and $\mu \in [0, 1]$ represents the viscous damping coefficient. The physical force $\mathbf{F}_i$ combines elastic repulsion and cadherin cohesion:
$$ \mathbf{F}_i = \sum_{j \neq i} \mathbf{F}_{\text{inter}}(i, j) $$
Let $r_{ij} = \|\mathbf{x}_j - \mathbf{x}_i\|$ represent the Euclidean distance between two cell centers, and $d_{ij} = R_i + R_j$ represent the sum of their radii. The interaction force vector $\mathbf{F}_{\text{inter}}(i, j)$ is computed along the unit normal direction $\hat{\mathbf{r}}_{ij}$:
$$ \mathbf{F}_{\text{inter}}(i, j) = \begin{cases}
-k_r (d_{ij} - r_{ij})\hat{\mathbf{r}}_{ij} & \text{if } r_{ij} < d_{ij} \quad (\text{Repulsion}) \\
k_a (r_{ij} - d_{ij})\hat{\mathbf{r}}_{ij} & \text{if } d_{ij} \le r_{ij} < 2.5 d_{ij} \quad (\text{Cohesion}) \\
0 & \text{if } r_{ij} \ge 2.5 d_{ij}
\end{cases} $$
Evaluating interaction forces across hundreds of independent cells would scale at $O(N^2)$ computational complexity, degrading rendering performance. To sustain high frame rates on consumer devices, we implement a localized spatial partitioning hash. The physical space is divided into discrete 2D spatial cells of width $w = 2.5 \cdot \max(R)$. Each cell's location $(x, y)$ is mapped to a structural 1D key:
$$ \text{hash}(x, y) = \left\lfloor \frac{x}{w} \right\rfloor + \left\lfloor \frac{y}{w} \right\rfloor \cdot W_{\text{grid}} $$
Collision checks are strictly restricted to the home cell and the surrounding contiguous 8 grid neighbors, capping collision checks to $O(N)$ computational complexity. The system features procedural real-time sonification via the browser Web Audio API, translating cellular state-changes dynamically into synthesized sound. Mitosis events synthesize high-frequency sine waves with an exponential sweep, while apoptotic necrotic conversions trigger low-frequency sawtooth sweeps, enhancing the multi-sensory educational output.
4. Future Directions
Future updates are focused on extending the platform's dimensionality and biological accuracy. Integrating a continuous 2D fluid dynamics solver utilizing the Navier-Stokes equations will allow us to simulate bioreactor fluid shear stresses and oxygen delivery profiles:
$$ \rho \left(\frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u}\right) = -\nabla p + \mu \nabla^2 \mathbf{u} + \mathbf{f} $$
This modeling addition will allow users to study dynamic media perfusion and how nutrient gradients impact necrosis. We are also designing complex tissue adhesion mechanics using Viscoelastic Kelvin-Voigt mechanical spring-dashpot models:
$$ \sigma(t) = E \varepsilon(t) + \eta \frac{d\varepsilon}{dt} $$
This approach will help capture the sticky mechanical properties of the surrounding extracellular matrix (ECM). Users will be able to model non-spherical structures, morphogenetic budding, and mechanical structural collapse.
Context-Aware Cross-Linking
If you wish to explore related scientific simulations, explore these tools: