Program Title: PolyMorph
CPC Library link to program files: TBA
Licensing provisions: BSD 3-clause
Programming language: C++11
Supplementary material: Figure fig:fig:1?
Journal reference of previous version: Comput.Phys.Commun. (2024) 109128, https://doi.org/10.1016/j.cpc.2024.109128
Does the new version supersede the previous version?: No
Nature of problem: In tissue development and disease, morphogenesis and cell fate determination depends on mechanical processes as well as chemical signaling. PolyMorph couples the Newtonian mechanics of deformable cells (including growth and proliferation) in 2D with a customizable set of reaction-advection-diffusion equations to simulate problems that require an integrated approach with chemical-mechanical interactions. Typical use cases include the patterning of epithelial tissues with chemical signals (e.g., morphogen gradients or the Turing mechanism), chemotaxis and cell migration, wound healing, diffusion- or nutrition-limited growth, regulatory network dynamics in a spatial cellular environment, and other problems in tissue self-organization. PolyMorph enables the numerical solution of such problems with bidirectional feedback between mechanics and chemistry, in large monolayer tissues and with an arbitrary number of interacting species.
Solution method: The off-lattice polygonal representation of cell boundaries in PolyHoop [1] is coupled to a lattice representation of diffusing chemical reactants. The reaction-advection-diffusion problem is solved with the finite difference method using the standard 5-point central difference stencil, and explicitly integrated in time. A scatter-gather approach inspired by the particle-in-cell method interpolates between the Lagrangian cell boundaries and the Eulerian finite-difference grid. The program is parallelized with OpenMP and does not use any external libraries.
Reasons for the new version: The original program PolyHoop solves the Newtonian mechanics of deformable particles, foams and cellular tissues in two dimensions. In the application realm of developmental and systems biology, however, various morphogenetic problems depend not only on cell mechanics, but also on chemical signaling, nutrient supply etc. The solution of such external transport problems was not part of PolyHoop, but is available in some related codes [2–5]. We have extended PolyHoop to create a new program named PolyMorph, introduced here. It enables the simulation of cellular tissue dynamics with bidirectional coupling to a multi-component diffusion system, which substantially widens the potential scope of application.
Summary of revisions: PolyMorph introduces a solver for \(n\) coupled reaction-advection-diffusion equations for concentrations \(c_i\), of the form \[\frac{\partial c_i}{\partial t}=D_i\left(\frac{\partial^2}{\partial x^2}+\alpha_i\frac{\partial^2}{\partial y^2}\right)c_i-\nabla\cdot(\vec{v}c_i)+R_i(p,c_1,...,c_n,k_1,...,k_m,t).\] The reaction terms \(R_i\) can be user-specified in the form of C++11 lambda functions, parameterized by the local polygon (cell) \(p\) and a set of \(m\) kinetic coefficients \(k_1,...,k_m\). Both the kinetic coefficients \(k_j(x,y)\) and the diffusivities \(D_i(x,y)\) (the latter assumed to be piece-wise constant) can access tissue information such as the local cell type through a map \((x,y)\to p\). This allows for cell-specific transport and reaction kinetics. At their birth, cells can draw their own random kinetic parameters to represent cell-to-cell variability [6,7]. In interstitial space or outside of the region occupied by the tissue, background coefficients can be specified. Note that the general form of \(R_i\) also allows for 0^th- or 1^st-order reactions, such as morphogen production or degradation. The coefficients \(\alpha_i\) control the degree of anisotropy in the diffusion of each species.
To simulate advection and dilution in moving, deforming and growing tissues, a discretized velocity field \(\vec{v}(x,y)\) is interpolated on the finite-difference grid from the underlying off-lattice motion of cellular boundaries. Inside cells, inverse distance weighting (IDW) of the cell vertex velocities is used. In the extracellular space, three options are available (Supplementary Fig. fig:fig:1?): i) IDW within a user-defined cutoff radius, ii) bilinear interpolation or iii) zero velocity. At the lattice borders, Dirichlet and Neumann boundary conditions can be specified.
For the coupling in opposite direction, all relevant cellular model parameters of PolyHoop, such as cortical tension, growth rates, etc., can be made dependent on a local readout of species concentrations \(c_i\) and their gradients \(\nabla c_i\), via user-specifiable C++11 lambda functions \(f(p,c_1,...,c_n,\nabla c_1,...,\nabla c_n,t)\). Contact parameters (the adhesion strength and the friction coefficient) can be defined based on the properties of both cells involved in the contact, allowing to model differential adhesion for instance. In addition, a chemotactic force vector and a cell type specifier can be defined in the above functional form.
Additional comments including restrictions and unusual features: PolyMorph inherits all physical features of its predecessor PolyHoop except cell fusion, but relaxes its design paradigm of strictly minimal code somewhat, in favor of more modularity to account for the increased level of complexity. It does therefore not supersede PolyHoop, but rather represents a spin-off program. The name PolyMorph conveys the shifted focus on morphogenetic problems involving a coupling to multi-component species transport.
Like PolyHoop, PolyMorph writes a series of VTK output files that can be viewed in ParaView. A structured grid file (.vts) containing user-specified lattice data is written for every animation frame. The user may further call the function Ensemble::write_OFF() to save a tissue in its current state in the common Object File Format (.off), for later use as a starting point of a different simulation.
Three limitations of the current implementation are that only diagonal elements in the diffusion coefficient matrix are supported, both the cell mechanics and the finite-difference solver operate with the same timestep, and both spatial directions use the same grid spacing.