|
LBIBCell
|
|
LBIBCell
|
Here, we explain the basic structure of an empty BioSolver. We recommend that, if you wish to implement a new custom BioSolver, that you start from this empty BioSolver and make sure that it is properly working. Once this is the case, the specific behaviour can be implemented. In the next sections, the implementation of various biological behaviour is discussed.
Let us now have a look at BioSolverEmpty.hpp. The header is surrounded by standard include guards:
Next, the BioBaseSolver.hpp is included, together with a few standard headers:
The BioSolvers are within the LbmLib::solver namespace:
Now it is time to declare our class, which inherits from BioBaseSolver and recursively uses itself as a template parameter:
We have to implement one single virtual public member function. This function is called at every timestep, and executes your custom biological behaviour. Since both a GeometryHandler and a ForceSolver object are passed, you have full access to the geometry, the fluid and CDE lattices, and the forces.
Then we have a little more boilerplate code: the BioBaseSolver has to be a friend class, the default constructor has to be private to prevent instantiation, and our BioSolver needs to have a unique name. Finally, the namespaces and include guards are closed.
Let us now discuss the implementation of our new class, which can be found in BioSolverEmpty.cpp. Be sides a few standard headers, we need to include our header, and the header to get access to the global simulation paramters (such as domain size, current time step):
For convenience, we define our locally used simlation parameters in an anonymous namespace (but this is a matter of taste). We recommend to define the frequency of applying this BioSolver:
Our private default constructor is empty:
Here we are: the function implementation which executes your custom biological operations. The first thing is to check whether it is time to execute the function (which is defined by FREQUENCY); if not, return without having done anything. In this case, the BioSolver simply prints a string to your terminal.
Finally, we have to give an unique name. In your application, you can add this BioSolver by simply knowing its name.
Finally, we need to get it run. This reqires only two simple steps. First, you need to include the two files (the header and source file) in your project. To be consistent, we want to add the files into the directories libs/LbmLib/include/solver/BioSolver/ and libs/LbmLib/src/solver/BioSolver/ (but any other place works as well, of course). Therefore, we have to add the following lines to the libs/LbmLib/CMakeLists.txt:
The last step is to add your BioSolver to your application. For instance, in tutorial_01.cpp, we would add the following line of code (you'll figure out the right location when reading Tutorial 1: Division, Differentiation and Signaling):
That's it, your new BioSolver works. In the following sections, we will not discuss the boilerplate code anymore, but only focus on the implementation of applyBioProcess(...).
For a general introduction on how to create a new BioSolver, please consult the section An Empty BioSolver. The implementation can be found in tutorial_01_BioSolverAreaRegulator.hpp and tutorial_01_BioSolverAreaRegulator.cpp.
First, we define some constants in an anonymous namespace which we will use later. We want to execute this solver at every time step. The target area of regulated cells shall be 200 units. The proportional constant KP defines how quickly the cell area is pulled towards the target area. Finally, MAXSOURCE defines an upper threshold for the mass source.
Now the implementation of the function applyBioProcess(...) is discussed. We specify two temporary maps. As keys, they store the domain identifiers of the cells, and the respective cell area (or mass source) will be stored as the value.
Now we compute the areas of the cell. This is actually implemented as a member function of the GeometryHandler object. The result is stored in our temporary areas map:
Now, we populate the target area map. This map (which is a private member, please see tutorial_01_BioSolverAreaRegulator.hpp) also stores the domain identifier flag of the cells as key (like the areas and sourcemap maps), and the target area as value. Here, all the cells get the same target area TARGETAREA. Cells of type 1 are later set to proliferate (this behaviour could already set here if desired).
Now the sources needed to maintain the desried target area are computed. The sourcemap tells us later what mass source a cell needs to converge towards its target area. This is simply done by computing the difference between the target and actual area, multiplied by fhe factor KP (note that more sophisticated controllers, e.g. PID controllers, could be easily implemented here):
Now we want to modify our mass source map a little bit. In particular, we do not want to exceed a certain maximal mass source (e.g. to limit the maximal growth rate). Here, we allow for negative mass sources, so cells also have the ability to shrink.
Let us loop the entire computational domain now in order to set the previously computed mass sources. Since we can do this in parallel, we can use openMP for the outermost loop:
For every lattice site, we now decide what mass source we apply. If the domain identifier is zero (that means we have surrounding or interstitial fluid), we do not add anything. For the cell type is 1, we add the maximal mass source MAXSOURCE, which leads to maximal cell growth and thus maximal proliferation. For the other cells, the area regulation kicks in: we use the previously computed mass sources.
For a general introduction on how to create a new BioSolver, please consult the section An Empty BioSolver. Here, only the implementation of the function applyBioProcess(...) is discussed. The implementation can be found in tutorial_01_BioSolverMembraneTension.hpp and tutorial_01_BioSolverMembraneTension.cpp.
Like already described above, we put our local constants into an anonymous namespace. We only apply this BioSolver every 10th time step; FORCE is the constant membrane force.
Now let's discuss applyBioProcess(...) member method. First, we get all the connections and store them in the tempconnectionmap:
To renew all the membrane forces, we first have to remove the old ones. The membrane forces have type 6 (please see apps/config/force.txt for a list of all available force types; you can also add your own force types if you wish to do so), so we remove all type 6 forces:
Finally, we loop all our connections and create new forces. Note that we have to create pairs of forces, since one force only acts on one geometry node, so we also have to create the mirrored force. The forces are added by creating a stringstream with the corresponding attributes (depending on the force type; please see apps/config/force.txt), which is added to the forcesolver:
For a general introduction on how to create a new BioSolver, please consult the section An Empty BioSolver. Here, only the implementation of the function applyBioProcess(...) is discussed. The implementation can be found in tutorial_01_BioSolverCellJunction.hpp and tutorial_01_BioSolverCellJunction.cpp.
As before, we put the local constants into an anonymous namespace. The solver is executed every 10th iteration. RADIUS is the radius within which a membrane point is looking for another memrbane point of another cell to create a cell-cell junction. If a cell-cell junction is formed, the spring constant will be K and the resting length L0.
Similar to the tutorial_01_BioSolverMembraneTension, we have to remove all the cell-cell junctions before new ones are created. Cell-cell junction forces are of type 7 (please see apps/config/force.txt for a list of all available force types; you can also add your own force types if you wish to do so), so we remove all type 7 forces:
The following lines of could could be simplified if you renounce parallel execution. Since openMP cannot handle std::map (and other non-random access containeres) directly, we have to first copy the content (the geometry nodes) into a std::vector:
Once this is done, we can loop all the geometry nodes in parallel:
Now we come to a central line of code: for node j, we are looking for other membrane nodes within RADIUS. For that, we use the getGeometryNodesWithinRadiusWithAvoidanceClosest member method, which returns the pointer to the closest geometry node - or nullptr if there is none.
If we have found a valid candidate, we gonna add a new junction. As for the tutorial_01_BioSolverMembraneTension, we actually have to add a pair of mirrored foreces, since one force only acts on one geometry node. The stringstream, which describes the force, is assembled according to the structure described in apps/config/force.txt. Finally, the new forces are added to the forcesolver.
For a general introduction on how to create a new BioSolver, please consult the section An Empty BioSolver. Here, only the implementation of the function applyBioProcess(...) is discussed. The implementation can be found in tutorial_01_BioSolverDifferentiation.hpp and tutorial_01_BioSolverDifferentiation.cpp.
As always, we put the local constant in an anonymous namespace. The THRESHOLD is a constant which is used for the differentiation. To be precise, for every cell, we integrate the signaling modecule concentration over the cell areas. Whenever this value will drop below that THRESHOLD, we will differentiate that cell from type 1 to type 2.
In a first step, we compute the spatial integral of our signaling molecule over the cell area. To do so, we use the method computeAccumulatedDomainConcentrations. It returns a map, where the keys are the domain identifiers of the cells, and the values the corresponding integrated concentrations:
Next, we create a map which takes the domain identifier as key, and a vector with all the connections of that cell. By that, we achieve a sorting and convenient access of the connections according to their cellular affiliation.
The following lines of code are only executed when the BioSolver is called the very first time. It is basically the inital condition for the cell types: the fluid always defaults to 0, and every initially existing cell defaults to cell type 1. Of course, this can be changed if you need something else. Also note that, if you load your geometry from a vtk file, the CellTypeTrackerMap is already automatically initialized, and the following lines of code are actually totally superfluous and will never be executed.
Now we can loop the map with the integrated concentrations, and decide whether we gonna differentiate that cell or not. In this case, if the value drops below THRESHOLD, we will differentiate to cell type 2:
The final step consists of copying the modified (differentiated) cell types to the lattice, i.e. the cell type flag of each lattice site is updated accordingly:
For a general introduction on how to create a new BioSolver, please consult the section An Empty BioSolver. Here, only the implementation of the function applyBioProcess(...) is discussed. The implementation can be found in tutorial_01_BioSolverCellDivision.hpp and tutorial_01_BioSolverCellDivision.cpp.
We set FREQUENCY to 100, i.e. the cells are only checked for cell division every 100th time step. The maximal cell size is set to 380, which means that a cell exceeding this value will be divided.
The first step is to update the internal map (celldefinition_) which stores the domain identifier as key, and all the connections of that cell as value:
Only when executed the first time, we have to initialize the local cell type map. This map stores the domain identifiers of the cells as key, and the corresponding cell type flag as value (initializing all the cells with type 1 is maybe not what you want in the general case, so you might adapt this to your needs). Also note that, if you load your geometry from a vtk file, the CellTypeTrackerMap is already automatically initialized, and the following lines of code are actually totally superfluous and will never be executed.
As previously for the tutorial_01_BioSolverAreaRegulator, we now compute the areas of the cells. The result is stored in a map, with the domain identifier of the cells as key, and the corresponding area as value. Here, we overwrite the area for domain identifier zero (which is the surrounding fluid) in order to prevent division of the surrounding fluid.
We now loop all the cells, have a look at their area, and decide whether they are big enough to be divided:
If yes, the cell is divided. The called function divideCell is relatively complex, and thus written in a separate function (which will be discussed below):
Because divideCell just altered the connectivity (and introduced a new cell), we have to update the local cell definition map (remember: domain identifier as key, all the corresponding connections as value):
We also have to add the new cell to the local cell type map (remember: domain identifier as key, corresponding cell type as value). Since the map is passed by reference, we can make changes here: (checkout getCellTypeTrackerMap), we can make changes
Now, the method cureLattice guarantees that the topological are consistent, i.e. the newly formed connections have to be integrated, connected to the lattice, boundary conditions have to be set up, etc.
Finally, the cell types (we just added a few new cells) have to be copied to the lattice (because each lattice site has a cell type flag):
So now, as promised, we discussed the member method divideCell. A new domain identifier is generated (simply increased the highest value by one), and the "affected" connections are stored in a temporary container. The affected connections are those which intersect with the cell division line. Ideally, only two affected connections are found; if less, something went really wrong; if more, the cell might be concave, and the cell division might lead to 3 daughter cells (we also want to avoid that). The cell division line itself can be defined in different ways. The software is shipped with two choices: random direction, or perpendicular to the longest axis of the cell. Here, we implemented the latter choice, discuss the corresponding method getTwoConnectionsLongestAxis in the end of this section. Of course, you can implement any other custom behaviour, depending on your biological hypothesis.
Now that we now which connections to cut, we detach the corresponding geometry nodes from those connections...
... and create two new connections which represent the cut. The attributes (boundary conditions, cell type, etc. are inherited from the affected (old) connections).
Now we can permanently remove the old connections:
There is one more thing: the one part of the old cell, which became the new cell, has to be updated since its connections still carry the old domain identifier flag. So we just choose one cell and update the domain identifier to newCellID:
So, how do we find the affected connections if we assume that the mother cell is cut perpendicular to the longest axis? This is implemented in the member method getTwoConnectionsLongestAxis which we will discuss now. Note that we make use of boost::geometry here; this is not really possible, but good to know in case you want to use some more fancy functions of boost::geometry. We just declare are few temporary variables which will be used later:
Now we loop all the connections of the cell which shall be divided, and search for those two points which have the largest distance (note that you could also conduct a PCA here in order to find principle axes):
Now that we found the longest axis, we compute the perpendicular line (the line which is used to cut the cells). This line is assumed to go through the midpoint of the longest axis:
Finally, we can search for the intersections between the cut line, and the connections of the cell. Basically, for each connection, we check whether both point of that connection are on the opposite side of the cut line. So if point 1 is on one side, and point 2 on the other side, we found an intersection. Or the other way round. If we haven't found extactly two connections, we'll throw an exception.
Whew! That BioSolver is quite complex. Luckily, if you want to implement other rules how to divide a cell, you can copy-past most of the code into your new BioSolver. Most probably, you'll only have to adapt a few lines of code.
For a general introduction on how to create a new BioSolver, please consult the section An Empty BioSolver. The implementation can be found in tutorial_01_BioSolverAreaRegulator.hpp and tutorial_01_BioSolverAreaRegulator.cpp.
First, we define some constants in an anonymous namespace which we will use later. We want to execute this solver at every time step. The MAXSOURCE defines the mass source strength. If the R^2L concentration locally exceeds the THRESHOLD, then MAXSOURCE is added.
Now the implementation of the function applyBioProcess(...) is discussed. Nothing fancy. We immediately return if the solver is not applied at this time step.
Let us loop the entire computational domain now. Since we can do this in parallel, we can use openMP for the outermost loop:
Since we only want to add mass for cells, we get rid of the surrounding and interstitial fluid. For every lattice site, we now get the local concentration of tutorial_02_CDESolverD2Q5_R and tutorial_02_CDESolverD2Q5_L and temporally store them.
If the complex concentration exceeds THRESHOLD, the source is set to be hight, and very small otherwise:
Finally, the source is assigned to the lattice site.
1.8.6 
