Add subcellular element model

[MILeach]

As part of the BBR grant, we should add the subcellular element model to the core chaste library.

Existing Code

This work builds on partially complete code initially developed in 2013.

The original copy of the code is available on trac A working copy is available in the Chaste/sem-project repository

Development Tasks

• [ ] Update the existing project code in the working copy to build against current develop
• [ ] Add new required functionality
• [ ] Add tests for new functionality
• [ ] Ensure subcellular element code meets required standards & is sufficiently tested
• [ ] Create tutorial & documentation

[AlexFletcher]

A brief overview of the code structure is provided in Section 6.2.3 of Dan Harvey's DPhil thesis.

It looks like Dan's approach is to create a class NodeBasedCellPopulationWithSemId inheriting from NodeBasedCellPopulation, where there is a 1-1 correspondence between nodes and cells, but not to actually use any cell functionality like division, growth, or death, so that cells aren't really thought of as 'actual' cells. Two NodeAttributes are then used to keep track of whether each node is actually part of an 'actual' cell's nucleus (mIsParticle) and to keep track of which 'actual' cell it belongs to (mIsRegion).

I guess this was done to limit the amount of new coding required, and because for the use cases considered by Dan Harvey and others, cell behaviours like division weren't necessary. I think a more general approach, allowing for cell division and cell death, will require a new cell population that maintains a correspondence between each cell and multiple nodes (as is the case already in VertexBasedCellPopulation and - once it is in the develop branch - ImmersedBoundaryCellPopulation).

With this new approach, I can see a few minor edits needed to the existing mechanics code, e.g.

• in SemForce::AddForceContribution(), to find which cell owns each node, instead of calling GetRegion() we'd want to get its global index then call GetCellUsingLocationIndex() on the cell population;
• in SemParameterScaler::ScaleSpringConstant(), instead of calling GetNumElements() with an semIndex, we'd instead want to get the number of nodes owned by a cell (elsewhere in the Chaste code, we'd do this by interrogating some kind of mesh element associated with the cell, so it depends on whether we'd want to stick to using NodesOnlyMesh in this new approach).

Then more effort would be required to implement things like cells dividing, growing/shrinking, dying, and interacting with any PDE state variables.

[AlexFletcher]

Dan Harvey's full DPhil thesis: https://ora.ox.ac.uk/objects/uuid:95f50f05-9cf5-4c58-9115-aff7aabdfd6f/files/m25393982233b14b3892ac994f014e96e

[AlexFletcher]

Original SEM paper: https://iopscience.iop.org/article/10.1088/1478-3975/5/1/015002/meta

[AlexFletcher]

Relevant recent paper by Dokmegang et al: https://www.biorxiv.org/content/10.1101/2020.06.08.140269v4

[fcooper8472]

The other SEM paper: https://doi.org/10.1088/1478-3975/8/4/045007

[fcooper8472]

Brief update, with a more thorough update coming later:

The current state of the PR #111 is that the code compiles, but simulations do not yet "work". The TestSemBasedSimulation compiles and runs, and produces some output interpretable by ParaView, but does yet do the right thing.

There is currently a large quantity of commented-out code which I am in the process of working through and re-adding.

[fcooper8472]

Update on finding the right place in the Chaste class hierarchy to put the SEM code

I've spent some time identifying where to put the SEM code, and discussed with @jmpf on 22-July. The overriding principle is that SEM is conceptually very similar to the overlapping spheres, but with groups of nodes forming individual cells.

Therefore, it should be possible to leverage the vast majority of existing functionality including, for as a priority, all the MPI parallelism.

---

Mesh

Currently we inherit from Abstract Mesh, but all the parallel functionality is in Nodes Only Mesh. Here, the "nodes only" part refers to the lack of (tetrahedral) elements. Leaving aside the absurdity of the name, a better option might be a NodesOnlyMeshWithElements, or something along those lines.

---

Cell population

Again, the node based cell population contains all infrastructure for parallel computation. It's slightly trickier to know how best to slot in the SEM based population into this hierarchy, because the Node Based population bakes in the assumption that there's a cell per node. It might be worth keeping this cell population separate as there is less machinery to duplicate in this case.

---

Forces

Fundamentally, the Generalised Linear Spring Force should be applicable directly to the SEM cell population. It is just fancy springs between nodes. But, the abstract two body force ONLY works with abstract centre based cell populations, which seems to be the incorrect abstraction.

We only inherit from the GeneralisedLinearSpringForce by bypassing the AddForceContribution method, getting around the safeguards.

We should modify the AbstractTwoBodyInterationForce to make it more general.

---

Recommendations

1. Mesh should inherit from NodesOnlyMesh
2. Cell population should PROBABLY stay where it is in the hierarchy...
3. Force hierarchy should be modified to allow a more general conception of spring-type forces