A combined experimental and computational framework to evaluate the behavior of therapeutic cells for peripheral nerve regeneration

Abstract Recent studies have explored the potential of tissue‐mimetic scaffolds in encouraging nerve regeneration. One of the major determinants of the regenerative success of cellular nerve repair constructs (NRCs) is the local microenvironment, particularly native low oxygen conditions which can affect implanted cell survival and functional performance. In vivo, cells reside in a range of environmental conditions due to the spatial gradients of nutrient concentrations that are established. Here we evaluate in vitro the differences in cellular behavior that such conditions induce, including key biological features such as oxygen metabolism, glucose consumption, cell death, and vascular endothelial growth factor secretion. Experimental measurements are used to devise and parameterize a mathematical model that describes the behavior of the cells. The proposed model effectively describes the interactions between cells and their microenvironment and could in the future be extended, allowing researchers to compare the behavior of different therapeutic cells. Such a combinatorial approach could be used to accelerate the clinical translation of NRCs by identifying which critical design features should be optimized when fabricating engineered nerve repair conduits.

being limited in their availability. Recently, research in PNI treatment has focused on nerve repair constructs (NRCs), that combine therapeutic cells and biomaterials. When implanted in the injury site, NRCs can provide mechanical support, guidance cues and a growthpermissive environment to modulate regeneration (Carriel et al., 2014;Hsu et al., 2013;Schuh et al., 2018). Research into the design of NRCs has so far focused predominantly on the choice of biomaterial, cell type and proregenerative cues. For NRCs that include a cellular component, key aspects that require optimization include spatial distribution of embedded cells and long-term nutrient supply.
The conditions in the local microenvironment, particularly the native low oxygen levels, are a major determinant of the regenerative success of cellular NRCs that are often overlooked. Under physiological conditions the characteristic penetration length of oxygen in tissue is considered to be around 100-200 μm, depending on the proximity to blood vessels, cell type, and density (Carrier et al., 2002;Rouwkema et al., 2008). PNIs cause acute damage and disruption to microvascular networks thereby obstructing tissue perfusion (Lim et al., 2015). The resulting local tissue hypoxia and absence of neovascularisation can affect oxygen diffusion and distribution within implanted NRCs. As a result, the core of the NRC, which often lies at a distance beyond the diffusion distance of oxygen away from the nerve stumps, may become severely hypoxic compared to the ends of the NRC. Besides the physical characteristics of the biomaterial used, oxygen consumption rates are also determined by the type of cells embedded in the construct (Cheema et al., 2012;Magliaro et al., 2019;McMurtrey, 2015). For instance, pluripotent stem cells tend to have low metabolic rates, while progenitor and differentiated cells have higher metabolic rates (Teslaa & Teitell, 2015). Additionally, the consumption rate of three-dimensional (3D) cultures appears to change based on the cell seeding density (Magliaro et al., 2019;Patzer II, 2004;Sielaff et al., 1997).
The formation of oxygen gradients has been found to correlate with gradients in viable cell density (Lewis et al., 2005;Radisic et al., 2006) and increased metabolic loads of the embedded cells (Carrier et al., 2002). Besides cell viability, oxygen availability is also linked to vascular regeneration. Cells in hypoxic environments often respond by the activation of multiple proangiogenic pathways and the upregulation of factors that encourage new vessel formation (Fong, 2008;Hashimoto & Shibasaki, 2015). Therefore, the presence of oxygen gradients within an NRC can lead to the formation of growth factor gradients, such as vascular endothelial growth factor (VEGF), which in turn can result in a more distinct directional chemotactic response of migrating endothelial cells. This can have further implications for PNIs as VEGF expression and neovascularisation have been found to induce axonal regrowth and Schwann cell proliferation and promote neural regeneration (Cattin et al., 2015;Donzelli et al., 2016;Hobson et al., 2000). Improving our understanding of the impact of the local microenvironment on implanted cells is therefore beneficial for informing NRC design. Nevertheless, most in vitro studies assessing the behavior and proregenerative capacity of NRCs, have been performed at standard laboratory incubator oxygen concentrations.
This condition does not represent the local in vivo endoneurial oxygen tension that studies in rat sciatic and human sural nerves report to be around 3%-7% (Newrick et al., 1986;Tuck et al., 1984;Zochodne et al., 1994). For this purpose, here, the effect of in vitro oxygen conditions on cell survival, VEGF release, as well as oxygen and glucose consumption, were measured. To simulate the sorts of cellular biomaterial used in translational nerve tissue engineering research, differentiated neural stem cells (CTX0E03) at a range of cell densities were embedded in thin, stabilized collagen constructs, and subsequently cultured at a physiologically relevant range of oxygen concentrations as well as standard cell culture conditions. CTX0E03 cells were selected as they are human clinical-grade therapeutic cells with demonstrated potential as an allogeneic "off the shelf" cell source for peripheral nerve repair (Kalladka et al., 2016;O'rourke et al., 2018;Rayner et al., 2021;Smith et al., 2012).
The in vitro cellular biomaterial model allows us to explore NRC performance in a highly controlled environment, providing an insight into the behavior of therapeutic cells in the critical first 24 h after implantation. However, while in vitro experiments are invaluable in furthering our understanding of cellular behavior and in improving NRC design, the vast number of possible scenarios, including differences in cell density and distribution, biomaterial permeability, and anisotropy, that need to be tested can be prohibitive. To accelerate the design process, mathematical modeling can be integrated into experimental work to create an efficient and robust multidisciplinary workflow and allow for continuous improvement through an iterative process (Coy et al., 2018). Such tools can also be used to extrapolate data from in vitro studies to an in vivo repair environment, thus refining future in vivo studies.
To this end, we derived a cell-solute model, which comprises a set of coupled partial differential equations describing the spatial and the temporal evolution of the CTX0E03 population and its local environment, including oxygen and glucose consumption and VEGF release. This allows us to assess the spatial gradients that will be established in repaired nerves and exploit this information in construct design (Coy et al., 2020(Coy et al., , 2021. This type of mathematical model has been widely used to enhance the design of engineered tissues and tissue culture bioreactors (Cochran et al., 2006;McMurtrey, 2015;Rutkowski & Heath, 2002), as they are computationally cost-effective, rely on a limited set of parameters, while still capturing most of the underlying biophysics. To calibrate this model, we first performed a sensitivity analysis of its outputs that allowed us to highlight the hierarchy between the different parameters. We then used the experimental observations to assign representative values to these parameters by matching the predictions of the mathematical model against the obtained data.

| Fabrication of stabilized cellular collagen gels
Differentiated CTX0E03 cells (dCTX0E03) were used to create stabilized cellular collagen scaffolds. These were used to mimic the conditions in engineered neural tissue constructs. All gels were prepared using 80% v/v type I rat tail collagen (2 mg/ml in 0.6% acetic acid; First Link) mixed with 10% v/v 10× minimum essential medium. The mixture was then neutralized using sodium hydroxide (NaOH) and 10% v/v cell suspension was added to give cellular collagen at a series of cell densities (0.5-1.5 × 10 6 cells/ml of gel).
Next, 240 μl of the cellular collagen mixture was added to individual wells of a 96 well plate and the gels were allowed to set at 37°C for 15 min. Using RAFT absorbers (Lonza Bioscience) the gels were stabilized using plastic compression for 15 min, a process whereby a biocompatible absorbent material is placed upon the gel and absorbs interstitial fluid to generate a dense, robust hydrogel (Brown et al., 2005). The resulting compressed gels were then immersed in culture medium and incubated at 37°C in a humidified incubator for 24 h under different oxygen concentrations, chosen to reflect the range of oxygen concentrations in which cells would reside in vivo.

| Low oxygen 3D cell culture and oxygen monitoring
A hypoxia workstation and incubator (HypoxyLab, Oxford Optronix) was used for experiments requiring low oxygen conditions. Cell culture medium for hypoxic experiments was conditioned to the target oxygen concentration (1%, 3%, or 7%), which encompasses the rage of endoneurium in vivo measurements, for 2 h before use.
Cellular collagen constructs were cultured in the HypoxyLab at the desired oxygen concentration for 24 h. In situ dissolved oxygen within the constructs was measured using the integrated OxyLite™ (Oxford Optronix) monitoring system. Fiber-optic oxygen probes (Oxford Optronix) were inserted into the middle of the constructs.
The sensor probes were set to continuously measure oxygen partial pressure (5 samplings/min). The results were recorded using Labview software (National Instruments). Results are presented as partial pressure values in mmHg (e.g. 7.6 mmHg corresponds to 1%).

| CellTiter-Glo assay
Metabolic activity was examined by measuring adenosine triphosphate (ATP) as an indicator and generating a luminescent readout, using the CellTiter-Glo ® 3D Assay (Promega). Based on the manufacturers protocol, a volume of reagent equal to that of the culture medium was added to each experimental well and following a 30-min incubation at room temperature, 200 μl from the assay solution were transferred to a microplate and luminescence was quantified on a plate reader (Flx800, BioTek). The metabolic activity of cells was determined in culture by measuring the intensity of luminescence signals after 24 h.

| Live/Dead assay
To assess cell viability, cultures were stained using Syto 21/ Propidium Iodide (PI) (Sigma Aldrich). Syto 21 is a green, fluorescent nucleic acid stain that exhibits bright, green fluorescence upon binding to nucleic acids in both live and dead cells. In comparison, PI, which exhibits red fluorescence, cannot permeate viable cells as it reaches the nucleus by passing through disordered areas of dead cell membrane. Thus, using both dyes allows for the simultaneous staining of viable and dead cells.
For the Syto 21/PI staining, the medium was removed from the gels, which were then washed three times with 200 μl of medium (37°C).
Subsequently, 200 μl of Syto 21/PI solution (1:1000 dilution) was added and the plates were incubated for 15 min at 37°C before removing the Syto-21/PI solution. The gels were then washed briefly with 200 μl of culture medium. Finally, an additional 200 μl of culture medium was added to each gel before image acquisition. Images were visualized using a confocal microscope (Zeiss-LSM710, Carl Zeiss) with ×20 phasecontrast water immersion objective. High-throughput quantification of cell viability from 3D image stacks by adapting a readily available ImageJ protocol was performed.

| VEGF release
The concentration of secreted vascular endothelial growth factor-A (VEGF-A) post 24 h incubation was determined by an enzyme-linked immunosorbent assay (ELISA). The cell medium supernatant from the gels was collected, stored at −20°C and later analyzed with a VEGF-A sandwich ELISA kit (human and rat VEGF-A kits, RayBiotech) according to the manufacturer's protocols.

| Experimental data analysis
Normality was determined using a Shapiro-Wilk test. Two-way statistical analysis of variance (ANOVA) was conducted, followed by Bonferroni's multiple comparison test.

| Cell-solute mathematical model
The cell-solute model is comprised of a set of continuous diffusionreaction equations that describe the interactions between oxygen (c), glucose (s) and VEGF (v) concentrations and the cell population (n), within the in vitro well setup. These variables were selected as they reflect the potential effect of the local microenvironment on the viability of seeded cells and the expression of VEGF, both of which are fundamental for nerve regeneration and vascular regeneration. We consider the gel and the medium above within the well geometry ( Figure 1) as continuous materials with effective uniform properties. In the following sections, we start by describing the model equations in the gel and culture medium, followed by initial and boundary conditions.
We consider the transport of oxygen in the gel to be driven by molecular diffusion, modeled using Fick's first law, and assume that cells metabolize oxygen following Michaelis-Menten kinetics, as is commonly used in the literature for conditions where oxygen is the limiting factor (Haselgrove et al., 1993;Huang et al., 2011;Magliaro et al., 2019;Zhong et al., 2018). (1) where c g represents the oxygen concentration in the gel, n the local cell density in the gel, D c g , the diffusion coefficient of oxygen in the gel, M c the maximum oxygen consumption rate by the cells, and c ̅ the oxygen concentration for which oxygen consumption by the cell is half its maximal value. Oxygen is considered to diffuse freely in the medium so that where c m represents the oxygen concentration in the medium and is the diffusion coefficient of oxygen in the medium.
Next, the equation that governs glucose concentration within the gel can be written as where s g represents the glucose concentration in the gel. Here, D s g , is the diffusion coefficient of glucose in the gel, M s the maximum glucose consumption rate by the cells, and s ̅ the glucose concentration for which glucose consumption by the cell is half its maximal value. Glucose consumption by cells is assumed to follow Michaelis-Menten kinetics (Aubert & Costalat, 2005;Dienel et al., 2017) modified with an additional term to capture anaerobic metabolism given the anticipated local oxygen This extra term enables the impact of anaerobic metabolism to be captured while only introducing one new parameter into the model. Similar to oxygen, we assume that glucose diffuses freely in the medium so that: where s m represents the glucose concentration in the medium and is the diffusion coefficient of oxygen in the medium.
Next, we consider the VEGF as an unstable molecule secreted by cells. Although vascular endothelial growth factors are a family of polypeptides, in this study we focus on modeling VEGF-A, which is considered the key mediator of angiogenesis (commonly referred to as VEGF). We describe the VEGF concentration in the gel by: where v g represents the VEGF concentration in the gel, D v g , is the diffusion coefficient of VEGF in the gel, K the VEGF degradation rate and G is the production rate of VEGF. Given that the production rate of VEGF by dCTX0E03 cells is not defined in the literature, we used the experiments presented in Figure 4 to define a production rate of VEGF that considers upregulation under low oxygen conditions. The final relationship describing the dependence of VEGF production on the underlying local environment is given by where α represents the baseline VEGF production rate and β represents the VEGF production rate depending on oxygen. Further c τ is the hypoxic threshold for VEGF production and n τ represents a crowding factor for the cells. VEGF is also assumed to diffuse freely in the medium: where v m represents the VEGF concentration in the medium and D v m , the diffusion coefficient of VEGF in the medium.
Finally, the viable cell density is determined by the balance of cell proliferation and death, along with cell migration. In collagen gels, however, cell migration is negligible on the short timescales considered here (Ardakani et al., 2014) and thus neglected.
In addition, CTX0E03 cells are conditionally immortalized and thus, do not proliferate in the absence of 4-OHT. We describe cell death as an increasing function of cell density, to represent competition for space, and a decreasing function of oxygen and glucose concentration, to represent competition for nutrients so that where δ c controls the oxygen-related death, δ s controls the glucoserelated death δ 0 encompasses all other interactions. The rationale behind the choice of oxygen and glucose-related deaths terms is exactly the same as that for anaerobic consumption (Equation 3).

| Well geometry and boundary conditions
The model consists of an axisymmetric 2D geometry (rotational symmetry along the vertical axis of the well), that represents the well of a 96-well plate and is composed of two domains: (i) the cell-seeded collagen gel at the base of the well, and (ii) the volume of culture medium above it

| Parameter values
Initially, the bounds of the parameters included in the model equations were informed based on literature values (see Table 2). For some parameters, literature was either scarce or conflicting and often from a range of cell types and culture conditions that do not represent the specific setup here. Thus, bounds for some parameters were chosen based on our own experimental observations and conditions.

| Sensitivity analysis
The model defined by Equations (1) where with μ* i being used to detect input parameters that have an overall influence on the output, and σ being used to detect input parameter involved in interactions or nonlinearities.

| Optimization
We seek to minimize the difference between the model predictions and experimental measurements by defining the underlying parameters. Given the size of the parameter set, we choose to perform a global, heuristic optimization using a particle swarm method (Kennedy & Eberhart, 1995). The particle swarm algorithm seeks to find an approximate solution to the equation: where θ is the vector of parameters, X { } the set or subset of species of interest (i.e., oxygen, glucose, VEGF, and cell density) and where J describes the cost function which when minimized corresponds to minimizing the average difference between experimental measurements and simulations for a given set of species. Thus Ω X is the size of the experimental set for species X, F θ σ ( , )

X I B
the vector containing the corresponding predicted values for a given vector of parameter value θ and vector of initial and boundary values σ IB , and Y σ ( ) X IB the vector corresponding to the corresponding experimental measurements. This approach has the advantages of avoiding possible local minima, considering the hierarchy between parameters, and imposing very few constraints on the regularity of the cost function itself.

| Numerical simulations
The model (Equations 1-9) was solved numerically using finite volume methods in Python 3.7. Given the rotational symmetry of a culture well, the model is solved in two-dimensional (radial and axial).
However, since the slope of the well geometry is small (<0.2%), variations in the radial direction can also be considered negligible compared to the axial ones, effectively rendering the model onedimensional. As for the axial direction, we devised a nonuniform twopart mesh corresponding to gel and medium with a change in mesh cell density at the interface between the two domains. This is done to allow a finer resolution in the gel where gradients are steeper, while still capturing the interface between gel and medium exactly.  Table 2 with 40 trajectories on a four-level grid.
Finally, the particle swarm optimization was performed using another open source library, PySwarm (Miranda, 2018), with three meta-parameters (two acceleration coefficients c1, c2 to control the individual and collective behavior, and one inertia coefficient w to control a history effect). The implementation was split in two separate steps (1) optimization of the oxygen, glucose and cell density related parameters (which are mutually coupled) (2) optimization of the VEGF related parameters (independent of the other species). For each step, we use 20 particles and 1250 samples leading to 25,000 simulations. For step 1, c1 = 2, c2 = 0.2, w = 0.6 whereas for step 2, c1 = 2, c2 = 0.2, w = 0.7 yielded the optimum results and met the appropriate convergence criteria (as defined in (Jiang et al., 2007)). During step (1)

| Viability and metabolic activity of dCTX0E03 cells under different oxygen conditions
The potential use of cellular NRCs for the treatment of PNI is dependent on the ability of encapsulated cells to remain viable and maintain their therapeutic effects. Therefore, the proportion of viable cells and their metabolic activity were evaluated. Figure 2a,b illustrate the survival of target cells under low oxygen conditions. Reduced oxygen availability seemed to cause impairment in dCTX0E03 survival and metabolic activity, with the effect being more pronounced at the highest seeding densities. There was up to 10% reduction in viability and up to 26% reduction in metabolic activity after incubation at 1% O 2 for 24 h, when compared with normoxic conditions. Differences in cellular responses can also be observed between the higher oxygen tensions, although they are not as pronounced. Finally, cell death appears to be higher at high cell seeding density conditions, possibly due to competition for available nutrients.
We also explored the spatial distribution of viable cells within the constructs (Figure 2c). Results indicate that areas of greater viable cell density occur at the top of the gels, correlating with highest oxygen concentrations at the air interface (and lowest at the well base which is furthest from the oxygen source) (Cheema et al., 2007). cells. Interestingly, the oxygen concentration appears to reincrease after 12 h for cellular constructs cultured at 7%. Between 0 and 12 h, F I G U R E 2 dCTX0E03 cell survival and metabolic activity in stabilized collagen gels exposed to different oxygen conditions for 24 h. (a) Cell viability was calculated using live/dead staining and analysis of obtained optical sections. Syto21 was used to label all cells and propidium iodide to label dead cells. (b) Metabolic activity was assessed using the 3D CellTiter-Glo assay. Data expressed as means ± SEM. Significance levels were *p < 0.033; **p < 0.002; and ***p < 0.001 compared with normal culture conditions (19%). (c) Spatial variability in the viability of dCTX0E03 cells in stabilized collagen gels exposed to different oxygen conditions for 24 h (60 × 10 6 cells/ml density after stabilization). Data expressed as means ± SEM (n = 4 independent repeats, three samples per condition).

| Oxygen consumption characteristics in 3D constructs
ELEFTHERIADOU ET AL. | 1987 dCTX0E03 cells cultured at 7% oxygen exhibit similar consumption characteristics as at other oxygen levels, namely a fast, initial decrease of the oxygen concentration, followed by much lower decrease rates. After 12 h a recovery is observed which could be attributed to a shift in the equilibrium between oxygen metabolism and supply. As a proportion of the embedded cells die, total oxygen consumption decreases which in turns leads to a reincrease in local oxygen levels.

| Functional analysis of dCTX0E03 cells under different oxygen conditions
Oxygen bioavailability is also directly linked to energy homeostasis.
Lower oxygen levels compromise the function of mitochondria in generating cellular energy currency, ATP, through oxidative phosphorylation, which is the most efficient way of producing ATP from glucose. This causes cells to rely on glycolytic ATP generation. Figure 4a demonstrates glucose consumption during the experiments. The glucose utilization rate was higher at low oxygen conditions. The decrease was also more pronounced at higher cell seeding densities.
Finally, as illustrated in Figure 4b, subjecting cells to physiological stress through oxygen deprivation stimulates and subsequently increases the expression of VEGF. VEGF release was affected by the local oxygen levels and cell seeding density, however, the relationship between them was not straightforward. Activation of VEGF expression by hypoxia-induced stress was more prominent at mild to severe hypoxia (1%-3%). For 1%-3% ambient oxygen concentrations, upregulation of VEGF release appears to reach maximum levels at n = 20 × 10 0 6 and n = 31 × 10 0 6 cells/ml, where cells were found to be more active. This trend was reversed for mildly hypoxic and normoxic conditions.

| Mathematical model
We derived a cell-solute model for a well geometry, which needed to be further parametrized for the specific cell type used. The sensitivity analysis enabled the prioritization of parameters that contribute most to variation in model predictions as summarized in Figure 5, where the values μ* (x-axis) and σ (y-axis) capture the impact of each model parameter on output predictions, and identify which parameters contribute to couplings or nonlinear effects, respectively. Based on previous literature the individual parameters can also be classified in terms of (non-) linearity, (non-) monotony based on their individual σ/ μ* ratio (Garcia Sanchez et al., 2014).
The oxygen concentration is mainly affected by the diffusion coefficient in the media, the oxygen concentration for which consumption is half-maximal, the maximum rate of oxygen consumption, the oxygenrelated death rate, and the baseline death rate (Figure 5a). The maximal rate of oxygen consumption has the strongest influence on model F I G U R E 3 Oxygen levels in the center of (a) acellular (n = 1) or dCTX0E03-seeded constructs (60 × 10 6 cells/ml density after stabilization) at (b) 1% oxygen, (c) 3% oxygen, (d) 7% oxygen. Time zero refers to the time point when the probe was positioned in the gel. Data expressed as means ± SEM (n = 3). s v g τ τ , as described in the Methods. Table 3 summarizes the final set of parameters found. against experimental observations appears to be worse for n = 60 × 10 0 6 cells/ml than for the other initial cell densities, but the model predictions closely follow the experimental data points.
Regardless of the initial cell seeding density, the simulated concentration of VEGF released into the medium after 24 h is also in good agreement with the corresponding experimental data. The poorest fit was found to be for n = 60 × 10 0 6 cells/ml, especially for 3% ambient oxygen concentration. In the case of oxygen consumption, the model qualitatively reproduces the general trend of the experimental data. For instance, in the case of 1% ambient oxygen the broad shape of the oxygen consumption curves matches that of the data, but the rates of decrease appear to be much quicker than the experimental values would suggest is realistic. This could indicate that the oxygen metabolism term requires further refinement in the future.

| DISCUSSION
This study explores the behavior in vitro of therapeutic cells under physiologically relevant oxygen conditions, one of the major determining factors that affect the performance of NRCs in vivo. With regard to oxygen, local supply after implantation is expected to be limited, especially during the first days when neovascularisation has not progressed. dCTX0E03 cells were found to be vulnerable to oxygen conditions they are likely to encounter in situ. However, the reduction of cell viability was not as significant as expected based on previous literature.
Extending the duration of the experiments could provide further insights regarding the low long-term survival upon implantation observed in previous studies (Smith et al., 2012;Stevanato et al., 2009). Moreover, this discrepancy could be associated with the Original values are divided by the initial cell seeding density. Data expressed as means ± SEM (n = 4 independent repeats, three samples per condition). Significance levels were *p < 0.033; **p < 0.002; and ***p < 0.001 compared to normal culture conditions (19%). Lee et al., 2016). We also found a correlation between increased glucose consumption and greater VEGF secretion, although this has not been investigated for neural stem cells before.
Next we developed a predictive, cell-type specific and computa-  (Table 2) for other cell types. One noteworthy exception is the maximal rate of glucose F I G U R E 5 Morris sensitivity analysis results based on final (a) oxygen, (b) VEGF, (c) glucose, (d) cell density values in the center of gel after 24 h. Each point represents the mean absolute value μ* (x-axis) and standard deviation σ (y-axis) of the elementary effect of each parameter; the first is used to identify which input parameters have an overall influence on the output (i.e., oxygen, VEGF, glucose, cell density) and the latter can help identify which input parameters are involved in interactions or nonlinearities.
consumption M ( ) s , which is almost 10 times higher than previously suggested (Gu et al., 2016;McMurtrey, 2015). However, the rate of glucose consumption by differentiated human neural stem cells, in general, is not widely characterized. Undoubtedly, the model described in this study involves a high degree of simplification of what are in reality complex biological phenomena. Nevertheless, the simulations appear to capture the cellular responses and related trends correctly. Some differences between the model outputs and experimental results were detected, with the largest ones being for n = 60 × 10 0 6 cells/ml at 3%, 7% O 2 .
To test whether these discrepancies were due to parameter estimation, we ran the PSO algorithm 10 times and confirmed that variability in the predicted parameter values was insufficient to account for the differences between measured and predicted VEGF concentrations (data not shown). This indicates that these differences were due to biological mechanisms that are not captured in the current governing equation set. For instance, the influence of VEGF concentration on the viable cell density was neglected here, even though it has been shown to influence the survival of neural stem cells under hypoxia. Moreover, from the two nutrients examined in Finally, another aspect that was ignored when modeling VEGF production and release is the presence of different isoforms. Cells are able to express different VEGF isoforms as part of their physiological processes (Ara et al., 2010;Cain et al., 2014). Still, including multiple species of the same molecule would have drastically increased the complexity of the model and the number of unknown parameters.
Each of the isoforms displays unique decay and diffusion characteristics, possibly due to differential collagen binding and proteolytic release (Vempati et al., 2011(Vempati et al., , 2014. Differential VEGF binding to collagen may be important during the generation of VEGF gradients within the construct and its release in the local microenvironment.
Therefore, including this mechanism in the model may improve its ability to predict the temporal and spatial VEGF distributions.
The overall quantitative framework that we developed by combining experimental and theoretical approaches can enable researchers to simulate a wide variety of different engineered tissue configurations and obtain robust predictions about the therapeutic effect of CTX0E03 cells embedded in NRCs. For instance, during the first critical hours upon implantation, therapeutic cells adapt to their environment by rapidly consuming oxygen. We could hypothesize that once the oxygen concentration reaches a value around the hypoxic threshold, the cells experience oxidative stress and produce VEGF that will later promote the migration of endothelial cells and neovascularization. This will in turn help perfuse the construct with oxygen and nutrients, supporting both the therapeutic cell population and the subsequent neurite outgrowth. Therefore, if we optimize the construct by identifying the design that yields the maximal viable cell density and most favorable VEGF gradients, we could potentially accelerate nerve regeneration. Moreover, if in the future a more comprehensive database of cell and material-type specific parameters is collated by repeating the in vitro experiments using different cell types, the mathematical model can be extended, allowing researchers to compare the behavior of different therapeutic cells under the same NRC configurations.

ACKNOWLEDGMENT
We are grateful for funding from a UCL Mechanical Engineering PhD Studentship (DE) as well as the EPSRC (EP/R004463/1).

CONFLICTS OF INTEREST
The authors declare no conflicts of interest.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are openly available.