# A 3D Agent-Based Model of Lung Fibrosis
## 2 Materials and Methods
### Extracellular Mediators  
As mentioned above, we implement an hybrid agent-based model, which is a model that encompasses both entities, called agents and fluids, that can diffuse through the simulation volume. More specifically, we include in our simulations extracellular mediators that affect the agents' environment. Mediators follow the general diffusion equation  
\[\frac{\partial C}{\partial t}-D_{\text{C}}\
abla^{2}C=S-\mu C \tag{1}\]  
Figure 2: Graphical representation of the simulation space as a dichotomic tree structure. where \(C\) is the concentration of the mediator, \(D_{C}\) is the isotropic diffusion coefficient, S is the source term (the cell secretion in our model), and \(\mu\) is the decay coefficient. The given coefficients (see Table 1) are specified at the beginning of a simulation for each substance together with the initial concentrations, which we set to be flat across the whole simulation space. In order to solve the partial differential equations (PDEs) in the form of Equation (1) for each of the substances, BioDynaMo implements multiple resolution schemes (we adopt the central difference method). To this end, both space and time are discretized: the diffusion grids (defined for each substance) split the simulation volume into evenly-spaced 3D boxes (whose size \(\Delta x=\Delta y=\Delta z\) can be defined by the user), while a constant global \(\Delta t\) defines the size of a time step. At each time step, the discretized concentration of each substance is updated within each box. In our model, we set the resolution (that is, the number of diffusion grid boxes along each dimension) to 10 for each substance so that (i) the Courant-Friedrichs-Lewy (CFL) stability condition \(\left(\frac{D_{c}\ast\Delta t}{\Delta x^{2}}\right)<\frac{1}{6}\) for the PDE solver is fulfilled (a necessary condition for the explicit version of the central finite difference method) and (ii) each box (with volume 1000 \(\upmu\)m\({}^{3}\)) can enclose an agent (with diameter 700 \(\upmu\)m). Moreover, we select the scheme \"DiffuseWithOpenEdge\" to make BioDynaMo run our simulations with open boundaries. In order to support interactions among different substances (which are needed by our model as it simulates ECM and TIMP depletion by MMP and MMP depletion by TIMP), we extended BioDynaMo by implementing this feature in the diffusion source code. Specifically, we edited  
* the function used to define substances at the beginning of a simulation so that both the depleting substance and the binding coefficient can be specified;
* the function that implements the central difference method by embedding the local depletion (i.e., \(\mu_{AB}\ast A_{i}\), where \(\mu_{AB}\) is the binding coefficient and \(A_{i}\) is the concentration of the depleting substance in the i-th box where the calculation is performed) into the decay term. To assess whether our custom BioDynaMo is able to properly simulate interactions among different substances, we develop a basic model involving MMP and TIMP. The mediators are secreted by stationary agents located at the center of the diffusion grid boxes, with one agent per box. Our model is then used to simulate the diffusion and interaction of the substances for a time frame of 2 days with different \(\Delta t\) values: 1 s, 0.1 s, 0.01 s, 0.001 s. As expected, the results presented in Figure 3 show a negative correlation between the size of the time step and the accuracy of the model. In fact, smaller \(\Delta t\) values lead to greater agreement between the steady-state reached by the model and the theoretical results (represented by the dashed red line).