Visualising cylindrical agents is not a simple task. However, since matrix fibres are typically much longer than their girth (for example, collagen fibres average 75 μm in length but only 2 μm in radius), we make the choice to visualise the cylinders as lines, at least for this release. We added two functions in the PhysiMeSS addon: fibre_agent_SVG and fibre_agent_legend. These functions are replacing PhysiCell’s default drawing function and instead draw ECM agents as lines in the outputted SVG plots and legend. In addition, once fibre cross-links have been calculated (after at least one mechanics timestep), the function fibre_agent_SVG colours the fibres according to the number of cross-links (more details in the section Determining fibre cross-links).

In Figure 2 we display the results of three simple 2D fibre arrangements, which can be reproduced using the .xml settings file Fibre_Initialisation/mymodel_initialisation.xml. In all three cases, we attempt to initialise 2,000 fibres of length 75 μm and radius 2 μm randomly (there is no .csv file prescribing their centres) across an 800 μm × 800 μm domain. In the left panel, fibres are placed isotropically (anisotropic_fibres is false) and so are given randomly assigned orientations. In the middle panel, fibres are placed anisotropically (anisotropic_fibres is true) with a prescribed orientation angle of 0.2 radians, so all fibres are aligned at an angle of 0.2 radians from the horizontal line. In the right panel, fibres are placed anisotropically with a prescribed orientation angle distributed around a mean of 0.2 radians with a standard deviation of 0.15 radians. We note that during the fibre setup, if a fibre is not contained wholly within the domain, it is disregarded. If the fibrous mesh is isotropic, we attempt to re-initialise such disregarded fibres, up to a maximum of ten times, by giving them a new orientation within the domain. If, after which time, the fibre still does not sit wholly within the domain, or the fibrous mesh is anisotropic (in which case changing the orientation will not be possible), it is permanently disregarded. This is why, in all three cases shown in Figure 2, the actual number of initialised fibres is less than 2,000. We note that the algorithm that disregards fibres not wholly contained within the domain leads to differences in fibre density at the boundary and may cause preferred fibre alignment with the boundary. For the cases we consider, we are interested in behaviours in the centre of the domain and far from the boundary; thus, this does not matter. Should the reader wish to study behaviours close to the boundary, they may wish to modify this algorithm or indeed use their own customised algorithm for fibre initialisation.

One immediate aspect to overcome in modelling cylindrical agents is how to determine neighbours. Finding neighbours is important in ensuring the efficiency and effectiveness of the code; for example, mechanical interactions only occur between neighbours. Currently, cellular agents in PhysiCell belong to a single voxel, and neighbours are found by checking for agents in the same or adjacent voxels. Fibres are typically much longer than a typical cell diameter (for example, collagen fibres are 75 μm in length; cells are 10 μm in diameter), so it is likely that they will extend across multiple voxels. Thus, to find neighbours of fibres, we need to correctly find and register each voxel a fibre passes through. To achieve this, we implemented custom methods to find neighbours in PhysiMeSS. Since these methods are common to both cells and fibres in PhysiMeSS, they are defined in the PhysiMeSS_Agent class. The method for finding and registering a fibre to its voxels is register_fibre_voxels, which is implemented differently for cells and fibres. We then build our own list of neighbours by first finding all the voxels that an agent, cell or fibre sits within, as well as all adjacent voxels. Neighbours are then determined to be any agents that are registered to any of these voxels via the method find_agent_neighbors. After finding neighbours, we de-register the fibre from all voxels bar the voxel containing its centre using the method deregister_fibre_voxels. Further details of these processes are provided in Appendix C..

A key aspect of the ECM is that fibres cross-link to form the matrix network. Within PhysiMeSS, we consider, at a basic level, that two fibres cross-link at any point that they touch or intersect. Within the PhysiMeSS class, we determine whether a pair of fibres cross-link using PhysiMeSS_Fibre’s method check_fibre_crosslinks. This is primarily a geometrical problem and is explained in detail in Appendix C.. It is only necessary to check cross-links for pairs of fibres that are in close proximity; thus, PhysiMeSS_Fibre’s method add_crosslinks determines cross-links between neighbouring fibres by calling the function check_fibre_crosslinks only if a pair of fibres are neighbours. Once fibre cross-links have been calculated, we colour the fibres according to the number of cross-links. This colouring follows a gradient of blue, with non-cross-linked fibres coloured light sky blue, fibres cross-linked once coloured steel blue, fibres cross-linked twice coloured blue and fibres cross-linked three or more times coloured dark blue. This is observed in Figure 3 where the initialised ECM fibres, as in Figure 2, are now coloured according to the number of cross-links. As an aside, and not unsurprisingly, we note an increase in the number of crosslinks (more dark blue fibres) for ECMs with more anisotropy.

A further geometrical problem is determining the nearest point on a fibre from any other given point, or more precisely, the displacement vector from a fibre to a point. This is calculated in PhysiMeSS_Fibre’s method nearest_point_on_fibre and is used within the check_fibre_crosslinks method. It is also used to determine the nearest point from a cell centre to a neighbouring fibre which we use to inform mechanical and chemical interactions between cells and fibres, see the section Fibre Specific Mechanics. Determining the nearest point on a fibre is discussed in more detail in Appendix C..

There are four possible interactions to consider, namely, cell–cell interactions, cell–fibre interactions, fibre–cell interactions and fibre–fibre interactions, thus four different potential functions are considered. Cell-cell interactions are handled by PhysiCell as per the function add_potentials. We do not go into detail here; more details can be found in [11].

Cell-fibre interactions are cell interactions with neighbouring fibres which affect the velocity of the cell, i.e. describe how a fibre affects a cell. These interactions are contained within the method add_potentials_from_fibre of the PhysiMeSS_Cell class. Firstly, we include simple repulsion and adhesion between cells and fibres as per between cells, copying the behaviour in add_potentials. This prevents cells from being able to go through fibres. Secondly, we introduce fibre-directed movement of cells as per [17] since cells are known to track along matrix fibres [18, 19]. We describe this briefly here, for more details please see [17]. We assume that cells move in response to a close neighbour fibre in two directions; an additional adhesive force, parallel to fibre orientation and an additional repulsive force orthogonal to the fibre (see, e.g., [20]) allowing the cell to track along the fibre.

The additional adhesive force takes the form It is directed along the normalised direction of the fibre, $\vec{\hat{f}}$, and depends on the normalised scalar product between $\vec{\hat{f}}$ and $\vec{v}$, the velocity of the cell, i.e., proportional to the vector projection of $\vec{v}$ onto $\vec{\hat{f}}$ (Figure 4). Moreover, the force depends on an adhesion coefficient, 𝛼fibre, and on a threshold velocity, vmax, which limits the pulling effect of fibres. We note that this equation implies the need for $|\vec{v}|\leq v_{\max }$. In the PhysiMeSS code, these parameters are called vel_adhesion and cell_velocity_max, respectively. The additional parameter s > 0 can be used to model additional effects that might increase (s < 1) or decrease (s > 1) the pulling effect. Following [17], we use s = 1.

The additional repulsion force is modelled as friction exerted by the fibre and takes the form It is directed parallel to the cell’s current velocity, v, and is affected by the component of cell velocity orthogonal to the fibre, i.e., proportional to the vector rejection of $\vec{v}$ onto $\vec{\hat{f}}$ (Figure 4). The strength of this force is dictated by 𝛽fibre, a cell–fibre friction coefficient, and the exponent r > 0 can be used to model nonlinear effects that increase (r < 1) or decrease (r > 1) the repulsion forces. Following, [17] we use r = 2. In the PhysiMeSS code, this parameter is called vel_contact. The additional cell–fibre interaction force is computed as the difference of these adhesion and repulsion terms, F = F∥ − F⊥.

Fibre-cell interactions are fibre interactions with neighbouring cells that affect the fibre. These interactions are contained within PhysiMeSS_Fibre’s method add_potentials_from_cell. Any fibre–cell interaction can only take place if a fibre has no more than one cross-link with another fibre; otherwise, we assume the fibre is fixed in place, tethered by its cross-links.

A single fibre with no cross-links can be pushed and/or rotated by a motile cell. Whether pushing or rotation are turned on can be controlled by the user using the Boolean flags fibre_pushing and fibre_rotation. Cell-fibre pushing provides repulsion to the fibre, giving it speed and moving its centre from its initial position. Cell-fibre rotation maintains the position of the fibre centre and modifies the coordinates of its extremities by changing the fibre orientation. The new fibre orientation is based on moment arm magnitude, impulse, and fibre length. The moment arm magnitude is determined using the Euclidean distance between the origin and the point of impact. The impulse is obtained by considering the friction of the fibre (modified by the user with parameter fibre_sticky) and the speed of the interacting cell. The calculated impulse and fibre length are then used to determine the angular velocity, representing the rate of change of the fibre’s orientation. Finally, the new orientation is obtained by rotating the old orientation vector using trigonometric functions. As a result, a cell impacting on a fibre end far from the centre will change the fibre’s orientation faster than impacting closer to the centre.

The mathematical equations that correspond to code implemented to reproduce the fibre rotation upon contact with a cell, are the following:

Moment arm magnitude:

Impulse:

Angular velocity:

New orientation:

An example of fibre rotation is available in Cell_Fibre_Mechanics/fibre_rotating.xml, and can be visualised in Figure 5. The parameters controlling the cell–fibre mechanics are defined as custom data and expressed in the cell definitions of fibre-interacting cells, in the .xml file and detailed in Table 2.

Fibre-fibre interactions are fibre interactions with neighbouring fibres. For this release, we do not consider that there are fibre–fibre interactions. However, the method add_potentials_from_fibre is incorporated into the PhysiMeSS_Fibre class for future development.

As discussed in the introduction to this paper, cellular remodelling of the ECM is a key process and an important interaction to consider between cells and matrix fibres. For this release of PhysiMeSS, we focused on matrix degradation. As well as degrading fibres, cells are responsible for creating new fibres and remodelling the ECM. The process of fibre generation by cells is not currently part of this release of PhysiMeSS but will be incorporated into a future release.

The degradation of the ECM plays a major role in tumour growth and invasion, for example. The key enzymes involved in this process are MMPs. MMPs can be found both within the ECM, at the cell membrane level (membrane-type MMPs), or being secreted by the cell itself. Different MMPs have different ECM-component targets, but their main function is to cleave and degrade structural ECM proteins such as elastin and collagen [21].

Fibre degradation is modelled in PhysiMeSS as a simple process of removing fibres from the domain under certain conditions. Fibre degradation can be turned on or off by the user using the Boolean flag fibre_degradation, a user parameter in the .xml file. During cell–fibre interactions as controlled by add_potentials_from_fibre, a fibre may be degraded if fibre degradation is turned on and further conditions are met, as detailed below. The fibre is then flagged for removal, which, for simplicity, happens instantly.

The first condition that must be met in order for a fibre to be degraded relates to whether a cell is hindered by the presence of a fibre. This may either be because as a cell tries to migrate through the domain, fibres act as obstacles in the cell’s path, or it may be because as cells proliferate, they require additional space for the growing mass of cells. To target the first of these problems, the user can define a parameter fibre_stuck_time that determines how long a cell is stuck in a certain place before it will try to degrade/remove a fibre in its path. To determine whether a cell is stuck, we check as part of the function physimess_update_cell_velocity whether the magnitude of a cell’s velocity is below a given threshold defined by the parameter fibre_stuck_threshold. In order that cells only remove fibres in their path, we consider the dot product of the cell’s motility vector with the vector connecting the cell centre to the nearest point on the fiber; if this dot product is greater than zero, a fibre is said to be in a cell’s path. A simple example of a cell migrating through the ECM and degrading cells along the way is available in Fibre_Degradation/mymodel_fibre_degradation.xml, and represented in Figure 6.

To target the second issue, whereby cells degrade fibres to make space for proliferation, we defined a new class in custom.h, called PhysiMeSS_Fibre_Custom_Degrade, which is derived from the PhysiMeSS_Fibre class. This new class implements in custom.cpp a custom degradation function, as described in Appendix D.. In this custom implementation, we check whether the pressure experienced by a cell is above the threshold fibre_pressure_threshold. In order to indicate the pressure experienced by a cell, we have introduced an additional user Boolean flag that, when set to true, colours the cells by the pressure they experience using color_cells_by_pressure. If either of the constraints outlined are satisfied and fibre degradation is turned on, the cell will degrade neighbouring fibres at the rate fibre_degradation_rate. The custom implementation is controlled by a Boolean parameter, fibre_custom_degradation. This example is available in Fibre_Degradation/mymodel_matrix_degradation.xml and represented in Figure 7, where we compare degradations with and without this custom, pressure-dependent degradation.

All the parameters describing the fibre degradation by cells are controlled by custom data expressed in the cell definitions of fibre-interacting cells, in .xml file and detailed in Table 3.