In-Silico Modeling of Tumor Spheroid Formation and Growth

Mathematical modeling has significant potential for understanding of biological models of cancer and to accelerate the progress in cross-disciplinary approaches of cancer treatment. In mathematical biology, solid tumor spheroids are often studied as preliminary in vitro models of avascular tumors. The size of spheroids and their cell number are easy to track, making them a simple in vitro model to investigate tumor behavior, quantitatively. The growth of solid tumors is comprised of three main stages: transient formation, monotonic growth and a plateau phase. The last two stages are extensively studied. However, the initial transient formation phase is typically missing from the literature. This stage is important in the early dynamics of growth, formation of clonal sub-populations, etc. In the current work, this transient formation is modeled by a reaction–diffusion partial differential equation (PDE) for cell concentration, coupled with an ordinary differential equation (ODE) for the spheroid radius. Analytical and numerical solutions of the coupled equations were obtained for the change in the radius of tumor spheroids over time. Human glioblastoma (hGB) cancer cells (U251 and U87) were spheroid cultured to validate the model prediction. Results of this study provide insight into the mechanism of development of solid tumors at their early stage of formation.


Introduction
Understanding tumors has been recognized as one the most challenging problems in biology and medicine. Their behavior involves complicated molecular biology and correspondingly complicated dynamics. Although tremendous effort has been devoted to developing therapeutic approaches for tumor suppression, there is still a significant need for new insights into complex aspects of tumors such as genetic instabilities, therapeutic resistance, inter-and intra-tumor heterogeneity, etc. Mathematical modeling is one of the most powerful approaches in predicting different aspects of tumor progression. It can provide quantitative prediction of biological processes and help interpret complicated physiological interactions in the tumor microenvironment. Provided collaboration with experimentalists, mathematical modeling can generate practical mechanistic solutions. Mathematical oncology, for example, promises to allow quantitative determination of effective individual therapies [1]. Therefore, the cross-disciplinary approach of combining mathematical and biological models of cancer has great potential to enhance cancer treatment.
Primary malignant tumors originate from small numbers of cells which become highly proliferative through mutations occurring in oncogenes, tumor suppressors or DNA repair genes [2]. Although the genetic mutations eventually lead to formation of large and complex vascular tumors, they all go through an avascular (devoid of blood vessels) mode at their early stage of growth [3]. Having insight into this stage of growth is valuable in understanding the behavior of tumors at subsequent phases. Tumor spheroids are predominantly used in in vitro models of avascular tumor growth [4][5][6][7][8]. They are spherical aggregations of cancer cells which are supplemented with a controlled amount of nutrient concentration. Spheroid size, cell numbers, and fractions of hypoxic, quiescent and proliferative cells are easy to track, which enables researchers to quantitatively investigate the effect of various parameters on the tumors' behavior [9][10][11].
In the continuum approach, the variables are assumed to be continuous, and described by means of PDEs that incorporate growth and kinetic interaction between components, and diffusion for transport phenomena [18][19][20]. For instance, Burton developed one of the early mathematical models to determine the radius of tumor spheroids by analysis of oxygen distribution inside the tumor [21]. The model could generate a growth curve, which was fitted well by a Gompertzian function. As another example of studies that examined the oxygen dynamics of spheroids, Grimes et al. demonstrated that the growth of spheroids can be obtained using two key parameters: rates of oxygen consumption and proliferation [22]. Greenspan established a model focusing on the growth-inhibiting effect of chemicals produced by disintegrated dead cells at the core of the tumor [23]. He concluded that the growth pattern strongly depends on a balance between the inhibitory effect of chemical components and the proliferation of living cells. Byrne and Chaplain introduced a model that provides the growth of tumor spheroids in response to local nutrient concentration [24]. They assumed that the tumor is an incompressible fluid in which the inside motion is induced by local proliferation of cells. They also incorporated cell-cell adhesion by the Gibbs-Thomson relation, which maintains the tumor's compactness. They further studied asymmetric perturbations to predict the modes of instability of radially symmetric growth. Further studies considered the effects of pH and the level of oxygen and nutrient as key parameters for the tumor growth rate [25,26]. Anisotropic growth of avascular tumor spheroids was also modeled in the context of continuum mechanics in references [27][28][29][30][31]. These studies modeled the tumor as a hyperelastic material and used multiplicative decomposition of the deformation gradient to investigate spatial distribution of stresses. They highlight the role of mechanical stress on the growth of tumors. Additional continuum-based studies of solid tumor growth are reviewed in references [32][33][34].
On the other hand, discrete approaches can capture the effects of cellular response, signaling pathways, inter-and intra-cellular interactions and tumor microenvironment on the growth of solid tumors [17]. For example, the dynamics of avascular tumor growth is presented using a lattice Monte Carlo model in [35]. The model predicts the growth of spheroids under various nutrient supply conditions. Hybrid models are another approach to modeling the growth of tumors that combine both continuum and discrete approaches to allow descriptions of macroscopic environmental variables (such as nutrient concentration) as well as discrete biological interactions. In the case of solid tumor growth, the clinical-size morphology of tumors can be studied in hybrid models while the cellular pathways and subcellular interactions are also involved [36]. For more details about discrete and hybrid approaches, the reader is encouraged to refer to the following papers [33,[37][38][39].
Development of biological patterns, such as animal coat markings and evolution of the enteric nervous system (ENS), is another type of process that often involves cellular proliferation and reaction/diffusion of materials [40][41][42]. To describe such processes, the evolution of a domain boundary is incorporated into the mathematical modeling [43][44][45][46]. This boundary evolution of patterns is very similar to that of tumors, i.e., growth, from a mathematical point of view. In ENS, for instance, the motion of neural precursor cells is studied by solving a PDE describing migrating cells [47][48][49]. Another example is the study done by Simpson to derive an exact solution for a linear reaction-diffusion PDE on a uniformly growing domain [50]. The method was verified by comparing the solution with a numerical approximation. Although they successfully obtained the exact solution for the density function in a growing domain, their solution is restricted to one-dimensional linear and exponential domain growth in Cartesian coordinates, whereas in tumor growth the evolution of a tumor's boundary is not prescribed and is obtained as the solution in spherical coordinates.
Despite the rich literature on mathematical modeling of solid tumor growth, the initial stage of tumor formation is either neglected or not comprehensively modeled. At this stage, the interplay between adhesive force, due to cell-cell interactions, and repulsive force, due to local pressure gradient, defines the dynamic formation of tumor spheroids, which further affects the subsequent growth phases. To better understand the involved mechanisms, there is still a need to conduct more collaborative research involving biological and mathematical modeling approaches.
Here, we focus on the early phase of tumor spheroid formation, which is mostly missing from previous studies. In this initial phase of spheroid formation, the cell-cell adhesion forces are dominant and drag cells toward each other to make a compact aggregation. This causes the spheroid to shrink. Once cells start to proliferate, the raised concentration of cells within the spheroid produces a pressure which compensates for the adhesion forces. At this point, the balance between forces stops spheroids from shrinking and leads to monotonic growth. This phase carries on until necrosis occurs where the competition between mitosis and necrosis defines the growth pattern. This early transient behavior of cells in tumor spheroid formation is modeled by a reaction-diffusion equation coupled with an ODE. The effect of adhesion forces between cells are incorporated into the system by prescribing a constant boundary condition. This fixed concentration, the so-called relaxed concentration, implicitly models an equilibrium condition at which adhesion forces and intrinsic pressure are in balance.
Both analytical and numerical solutions were obtained for the change in the radius of a tumor spheroid at the early stage of formation. The theoretical model was validated against the formation and growth of tumor spheroids generated from glioblastoma cell lines.

Model Formulation
A tumor spheroid is considered to be a system of particles (cells) with continuous change in their properties such as concentration, velocity, etc. To start modeling this system, the material time derivative of the tumor spheroid's mass is written as a balance equation for the continuum concentration of cells within the system: where m is the mass of the tumor, C(x, t) is local concentration of cells, x is the position vector, Π is a volumetric source of mass, dV is a volume element and ∂ ∂t V C(x, t)dV is a material time derivative when particle locations are held fixed. Using the Reynolds transport theorem, one can write where v(x, t) is the velocity of particles. Followed by localization, Equation (1) gives: A constitutive statement analogous to Darcy's law can be proposed for the flux of cells as is the pressure inside the tumor and k is the proportionality of cell flux to the gradient of pressure, noting that cells move from higher pressure to lower pressure regions. Here we propose that this pressure is linearly proportional to the concentration of cells. Hence, Equation (3) yields: where K denotes ability of cells to respond to the gradient of concentration. This property depends on cell diffusivity (D) in response to the gradient of the concentration and adhesion forces between cells ( f ), i.e., K =K(D, f ). Cell adhesion reduces the ability of cells to move in response to the gradient of concentration. Therefore, to simplify the model we propose a linear relation, K = D, where 0 < < 1. The volumetric growth depends on the type of cells, concentrations of nutrients, growth factors and environmental cues. Here we propose a simple linear relation for volumetric growth, where η is the rate of cell proliferation. This assumption is correct for early growth of tumor spheroids since they are usually small, such that all cells can receive enough oxygen and nutrient and no hypoxia occurs in the tumor. Under this assumption, Equation (4) can be rewritten as: Spherical coordinates, which are well suited to represent the spherical shape of tumor spheroids, can be used to advantage. Since the early stage of growth is radially symmetric, our analysis focuses on radially symmetric solutions. Equation (5) in spherical coordinates become: subject to the following boundary and initial conditions: where R(t) is tumor radius, C 0 is the imposed boundary condition on concentration of cells, and C i is the initial cell concentration, assumed spatially constant. Denoting the free surface boundary by , the change in radius can be derived by integrating cell motion on the entire volume of the tumor as (see Appendix A.1) in which R 0 is the initial radius. The proliferation inside the tumor increases local concentration of cells and generates a pressure gradient, and consequently a concentration gradient. Equation (5) indicates that cells move from higher concentration regions to lower concentration regions to create a uniform concentration (relaxed concentration) everywhere so that cells do not feel any extra pressure. For this pressure increment to be stabilized by the adhesion forces between cells, the volume of the tumor must increase to reach the relaxed concentration. To model this equilibrium, we introduce boundary condition (7b) which imposes a constant relaxed concentration on the boundary of the tumor spheroid. Solutions to Equations (6) and (8) with corresponding initial and boundary conditions in (7) give the distribution of cells and change in radius of the tumor spheroids over time.

Analytical Solution
Equation (6) is a linear concentration-dependent reaction-diffusion equation with mixed boundary conditions. The reaction-diffusion form with constant source term is: is a solution of (9), the following is a solution of (6) (see Appendix A.2): If a solution to (9) satisfies the initial and boundary conditions (7), then so does (10). Therefore, the first step is to solve Equation (9) subject to initial and boundary conditions (7). Using the variable change C 1 (r, t) = U(r, t) + C 0 , Equation (9) becomes a standard homogeneous PDE with zero boundary conditions and constant initial condition, as follows: where R 1 (t) is the radius of tumor spheroids without proliferation. The following solution is obtained for The solution to the full reaction-diffusion Equation (5), is obtained as (see Appendix A.4):

Model Simplification
By introducing a new variable,R = nπ The timescale of proliferation is small enough compared to that of cell motility such that at each instant the radius of the diffusion-only model, R 1 (t), is close to the concentration-dependent reaction-diffusion one, R(t). In the limit of separation of time scales this approximation becomes exact. Therefore, it can be assumed thatR(t) ≈ nπ. Substitution of Equation (14) into Equation (13) simplifies the rate of change of radius to which, after simplifications, has the following solution (see Appendix A.5)

Numerical Solution
In this section, numerical solutions of the model expressed in Equations (6) and (8) are presented. Equation (6) is a reaction-diffusion equation with mixed boundary conditions coupled with the ODE in Equation (8). To solve this system of equations, temporal and spatial discretizations are required. The boundary of the tumor spheroid is moving in time, which requires a new spatial discretization at each time-step. To keep the number of nodes constant, the position of each node must be able to move in time. To this end, a mapping , where τ = t max and R(t) is the moving boundary of the tumor. Using this change of variables, the moving domain of the solution is mapped to a new domain which always stays between zero and one. Equations (6) and (8) can be rewritten in the new variables as where λ = τK R(t ) 2 , subject to the following boundary and initial conditions: Here to solve this system of equations, the Crank-Nicolson (CN) finite difference scheme was employed [51] (see Appendix A.6).

Model Analysis
The term α in Equation (16) contains the influence of both cell motility in response to the gradient of concentration, K, and the relaxed concentration of cells, C 0 . The higher the absolute difference, C 0 − C i , the more shrinkage is expected. This qualitative effect holds for K as well, i.e., the higher the motility, the faster cells respond to the gradient of concentration. To show this effect quantitatively, the value of the normalized radius, , is plotted in Figure 1a for C i C 0 and different values of K * = K K 0 holding other parameters fixed. Please note that we selected t max = 210 h to be consistent with the experimental results in the next section. Additionally, we take K 0 to be 10 −10 cm 2 ·s −1 [52]. As can be seen, the shrinkage of the tumor is faster for cells with higher motility (K). The tumor spheroid decays further until the concentration of cells reaches the relaxed concentration where the diffusivity of cells and adhesion forces are in balance (minimum tumor radius). At this point, the proliferation continues to elevate cell concentration and breaks the balance. To reach a new balance, the tumor spheroid increases its radius to reduce the local concentration which finally leads to monotonic growth. Rate of proliferation for highly proliferative cancer cells is normally in the order of η 0 = 10 −2 h −1 . This value was used as a reference number to non-dimensionalize the proliferation rates in our analysis. The effect of cell motility on formation of spheroids is illustrated in Figure 1a in which parameters are set as η = 1.8 × η 0 , C 0 C i = 1.5 and R 0 = 0.01 cm. The higher K in the figure corresponds to lower minimum radius and faster shrinkage. This result also shows that tumor spheroids with higher K grow faster since cells can rapidly respond to local proliferation and reach balance by moving the boundary of the tumor. Unlike many types of mammalian cells which have an intrinsic cell program that restricts their proliferation, most cancer cells are highly proliferative [2]. When cells proliferate, the local concentration increases, and the generated pressure moves cells away. The formation of a tumor spheroid is faster if cells have a high proliferation rate. To illustrate this effect, the formation of a tumor spheroid is depicted in Figure 1b for different values of η * = η η 0 , holding other parameters fixed, i.e., K * = 1, C 0 C i = 1.5 and R 0 = 0.01 cm. As can be seen, a tumor spheroid with a higher proliferation rate assembles faster.  Figure 1. Formation of tumor spheroids with different values of (a) cell motility, K * 1 = 1, K * 2 = 2 and K * 3 = 3 (with η * = 1.8, C 0 C i = 1.5 and R 0 = 0.01 cm), (b) proliferation rate, η * 1 = 1.6, η * 2 = 1.8 and η * 3 = 2 (with K * = 1, C 0 C i = 1.5 and R 0 = 0.01 cm). A tumor spheroid with higher cell motility grows faster since cells can rapidly respond to local proliferation and reach balance by moving the boundary of the tumor. Additionally, a tumor spheroid with a higher proliferation rate assembles faster.
To compare the analytical and numerical solutions, the formation of the tumor spheroids was obtained for two sets of parameters, (i) η * = 1.6, C 0 C i = 1.5, K * = 1, R 0 = 0.01 cm, and (ii) η * = 2, C 0 C i = 2, K * = 1, R 0 = 0.02 cm, using the Crank-Nicolson scheme outlined in Appendix A.6. Results are compared with the analytical solution in Figure 2.

Analytical Analytical
Numerical ( high proliferation rate-set1) Numerical ( low proliferation rate-set2 ) Figure 2. Formation of tumor spheroid obtained by analytical solution and numerical prediction, using two sets of parameters; set 1: η * = 1.6, C 0 C i = 1.5, K * = 1, R 0 = 0.01 cm, and set 2: η * = 2, In the contraction phase the analytical and numerical solutions reasonably match. The analytical solution loses accuracy once growth becomes dominant. This is a result of the simplification we made in Equation (14).
As can be seen in the figure, analytical and numerical solutions match very well in the contraction phase. The analytical solution loses accuracy once growth becomes dominant. This is a result of the simplification we made in Equation (14). Both solutions predict that the tumor spheroid in parameter set 2, which has higher proliferation rate, experiences faster contraction and faster growth.

Model Validation
In this section, the theoretical model is validated against the formation of in vitro solid tumor spheroids generated from glioma cell lines (U251 and U87 hGB cells), as the most lethal type of intracranial tumors. Reproducibility, ease of assembly and ability to provide high-throughput screening make them a promising candidate for in vitro threedimensional tumor models. Figure 3 shows the rates of proliferation of the two cell lines, η U251 = 0.037 ± 0.004 h −1 and η U87 = 0.026 ± 0.003 h −1 , which are in the range of data reported in [53], i.e., η U251 = 0.038 h −1 and η U87 = 0.033 h −1 . Results of U251 spheroid formation over 210 h are shown in Figure 4. During the formation phase, intercellular interactions generate adhesion forces which pull cells together and increase the concentration of cells within the spheroid. The size of the spheroid reduces since the proliferation is not yet dominant. The tumor spheroid shrinks until the concentration reaches the relaxed concentration (C 0 ). At this minimum radius, the driving forces are in balance, i.e., adhesion forces and forces due to high concentration of cells within the spheroid. This balance breaks once the proliferation of cells becomes dominant, elevating the local concentration above the relaxed concentration. To remove the produced force, the boundary of the tumor spheroid moves to increase the volume. This volume increment reduces the concentration of cells and equilibrates the forces inside the tumor spheroid.
For spheroids which did not have a full spherical shape, the average of the largest and smallest diameters was considered to be the spheroid diameter. Figure 5 shows the size of the tumor spheroids over time compared with analytical and numerical solutions. As shown in the figure, the mathematical model provides a reasonable prediction of the formation of tumor spheroids and the minimum diameter. The model predictions were able to follow the trend of formation and growth until approximately 160∼180 h, after which the tumor spheroids lost their homogeneity in terms of cell viability level. It is evident that in big spheroids, cells close to the core become hypoxic and change their metabolism. This can reduce the accuracy of the model.   The model is able to predict the formation of tumor spheroids and the minimum diameter, but loses accuracy after approximately 160∼180 h. This divergence from experimental results is started when the tumor spheroids lose their homogeneity due to hypoxia and/or necrosis initiation.

Proliferation Rate
HGB cells were cultured in Dulbecco's Modification of Eagle's Medium (DMEM) supplemented with 10% (v/v) Fetal Bovine Serum (FBS) and 1% (v/v) Penicillin/Streptomycin, and incubated at 37°C in 5% CO 2 . The number of cells, N 1 and N 2 , was counted using a Trypan blue assay after 24 h and 48 h, respectively (fresh media was added after 24 h). Rates of proliferation (1/h) were calculated, shown in Figure 3, for both cell lines as

Spheroid Culture
For culturing spheroids, cells cultured in Section 4.1, were dissociated with Gibco TM Trypsin-EDTA (0.5%) and were centrifuged at 300× g for 5 min. After removing the supernatant and suspending the cell pellet in 1 mL of medium, the number of cells was counted using a Trypan blue assay. Afterwards, self-filling micro-well arrays (SFMAs) were used to produce uniform tumor spheroids [54]. The desired concentration of cells was loaded dropwise through guiding channels of SFMAs and were gently seeded into the wells. The microwells were kept in an incubator and imaging started 5 h after seeding to let the cells fully settle in the wells. Cells were supplemented with fresh medium every 24 h to maintain the concentration of nutrients. The formation of spheroids was imaged using optical microscopy (Axio Observer, ZEISS, Oberkochen, Germany) over 210 h. The size of spheroids was measured using ImageJ [55].

Conclusions
In this work, we have presented an analytical solution for the formation of solid tumor spheroids. The process of tumor spheroid formation includes a preliminary contraction phase where adhesion forces densify cell aggregation. This phase proceeds until the cell concentration reaches a threshold, the so-called "relaxed concentration" at equilibrium. Afterwards, cell proliferation raises concentration and produces pressure which breaks the equilibrium. The tumor spheroid evolves in size to compensate for the generated pressure. This transient phase in formation and growth of tumor spheroids was mathematically modeled using a system of coupled PDE and ODE with appropriate boundary and initial conditions. To validate model predictions, human glioblastoma cancer cell lines were spheroid cultured and their size was imaged over 210 h. Results showed that although the model loses accuracy after approximately 160∼180 h, it can nevertheless provide reliable prediction of the size of the spheroids before they become inhomogeneous.It should be noted that this study is limited to the modeling of solid tumor formation with no access to environmental stimuli such as stroma, immune cells, etc. However, our approach has the potential to include the inhibitory effect of drugs using an additional reaction-diffusion equation. The effect of a drug on the tumor development may be expected either to simply extend the contraction phase, or to cause a monotonic contraction after the expansion, in either case due to the apoptotic effect of the drug. In this section, we derive an expression to relate the gradient of concentration at the boundary of a tumor to the time rate of change of its radius.
The rate of change of volume is where n is normal to the surface, i.e., n =ê r . Additionally, ∇C(r, t) = ∂C(r, t) ∂rê r , and . The surface element on a sphere of radius R is dΩ = R 2 sin(φ)dθdφ, where θ and φ are azimuthal and polar angles, respectively. Hence, Equation (A1) gives: which further simplifies to: . (A3)

Appendix A.2. Proof of Proposition 1
In this section, the proof of Proposition 1, which was used as an intermediate step to solve Equation (6), is provided.
The RHS of (6) is Hence Equation (6) is satisfied.
Here, we solve Equation (11a) and derive an approximate expression for R 1 (t).
Using separation of variables, U(r, t) = P(r) T(t), one can derive the solution of (11a) as where D n = 2R(0) (C 0 − C i ) (−1) n nπ and λ n = nπ R 1 (t) . Hence, Substitution of (A6) into (A3) gives The infinite series in Equation (A7) is convergent and has upper and lower bounds as follows: Using the identity e −(λ 2 Here, we observe that the upper and lower bounds in Equation (A9) are very tight, at least for system parameters in a reasonable range. Therefore, considering either bound as an approximate solution for dR 1 (t) dt is acceptable. To show this, the solution of R 1 (t) using the upper bound, R 1 (t) ub , and lower bound, R 1 (t) lb , are found by a finite difference method and the corresponding relative error, , is shown in Figure A1.
System parameters are adopted in the range of biological properties of tumor spheroids as K = 10 −10 cm 2 /s, R 0 = 0.01 cm. We also assume that the initial phase of growth does not take longer than a few days, t max =240 h, and As can be seen in the figure, the error is very small, indicating that both bounds are acceptable in this range. To reduce the complexity in the analytical solution procedure, we used the upper bound as the solution of dR 1 (t) dt , for which the analytical solution is available as which further simplifies to: It should be noted that there are two radii in this solution, i.e., R 1 (t) and R(t). The former is the radius of the reaction-diffusion form with constant source, Equation (9), and the latter is the radius of the concentration-dependent reaction-diffusion form, Equation (5). By substitution of (A12) into (A3), the first order differential equation for R(t) is obtained as: where F(R, n) = λ n R(t) cos(λ n R(t)) − sin(λ n R(t)) R(t) 2 .

Appendix A.5. Model Simplification
Here, steps to simplify Equation (15) are provided. Equation (15) gives the rate of change of radius as Using the same identity used in (A9) and the argument proposed in deriving Equation (A10), Equation (A14) simplifies to where α = f i √ πK 2 . It should be noted that C i K 1, (1 + C i K ≈ 1). Hence, Equation (A15) simplifies approximately to the following nonlinear integro-differential equation: Please note that R(0) is the initial radius, and R 1 (t) and R(t) are radii for the constantsource reaction-diffusion equation and the concentration-dependent reaction-diffusion equation, respectively. It can be shown that the term R(0)R 1 (t) R 2 (t) fluctuates around 1 during the shrinkage and early growth phase as follows: R(0)R 1 (t) R 2 (t) < 1, in the initial growth phase .
Thus, we write R(0)R 1 (t) R 2 (t) = 1 + ε(t), where ε(t) denotes the variation around 1 and depends on the intrinsic properties of the system, such as proliferation rate (η), cell diffusivity (D), etc. Although the form of function ε(t) cannot be specified precisely, we make the approximation ε(t) ≈ 0 to simplify Equation (A16). This assumption is not accurate, especially for the growth phase. However, the comparison between the analytical solution and both numerical and experimental results in Sections 3.1 and 3.2 confirms that this approximation provides an acceptable prediction on the formation of spheroids in the range of interest. Using the identity t 0 e ητ √ τ dτ = π η erfi( ηt), Equation (A16) then yields: Please note that if C i = C 0 then α = 0, and the size of the tumor does not change, but experimentally, the initial seeding concentration, C i , is always less than C 0 , since a compact tumor mass has not yet formed, so this can be considered physically implausible. Using the fact that the solution to (A18) can be obtained as Appendix A.6. Numerical Method The Crank-Nicolson scheme is implicit, unconditionally stable, and gives secondorder convergence. Considering temporal (i) and spatial (j) discretizations, the following approximations to temporal and spatial derivatives transform the differential Equations (6) and (8) into an algebraic system of equations.
∂C(i, j) ∂r Substitution of these approximations into (17) gives a system of linear algebraic equations in the form of Ax = b, where A is square matrix of coefficients and vector x contains unknown concentrations at each node. The CN scheme is implicit in both time and space. Based on this method, the (k + 1) th iteration of unknown concentrations, C (k+1) , is defined as: where A lower and A upper are respectively a lower triangular and strictly upper triangular decomposition of matrix A. As the iterations proceed, the approximations converge until the error reaches a defined tolerance. MATLAB [56] was used to perform iterations to obtain unknown concentrations at each node. Please note that we used explicit value of R(t) in Equation (17a), i.e., R i (t), such that we could solve for the concentrations first and plug them into Equation (17b) to obtain R i+1 (t).