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

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.

INTRODUCTION

Peripheral nerve injuries (PNIs) are associated with high socioeconomic and personal costs; the mean patient age is ~30 years so they can greatly impact lifetime health and productivity (Lad et al., 2010). Although nerves exhibit some regenerative capacity, the degree of reinnervation and subsequent recovery is dependent on many factors. In the case of incomplete regeneration, patients might not achieve meaningful functional recovery and the resulting disability may be highly debilitating with long‐term effects on them and their families (Grinsell & Keating, 2014; Panagopoulos et al., 2017).

The current “gold standard” of treatment for large gap PNIs is an autograft, where a sensory nerve is harvested and sutured to bridge the gap. However, autografts exhibit success rates that are far from ideal (Yang et al., 2011), as well as causing donor site morbidity and 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 growth‐permissive 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, 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.

METHODS

Unless otherwise stated all cell culture materials were purchased from Sigma Aldrich or Thermo Fisher Scientific.

Culture and differentiation of CTX0E03 cells

Human neural stem cells (CTX0E03, level P28‐P31, ReNeuron Ltd) were cultured in Dulbecco's modified Eagle's medium:F12 medium supplemented with human albumin (0.03%; Nova Biologics); Glutamax (2 mM); human transferrin (5 μg/ml), putrescine dihydrochloride (16.2 μg/ml), human insulin (5 μg/ml; Sigma), progesterone (60 ng/ml; Sigma), sodium selenite (40 ng/ml), epidermal growth factor (20 ng/ml), basic fibroblast growth factor (10 ng/ml; Invitrogen), and 4‐hydroxytamoxifen 4‐OHT (100 nM) in 175 cm2 laminin‐coated (10 µg/ml; Amsbio) flasks. Following expansion, CTX0E03 cells were subsequently differentiated for 1 week by removal of growth factors and 4‐OHT.

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 × 106 cells/ml of gel). These cell seeding densities were based upon the range used within NRCs (Coy et al., 2020; Georgiou et al., 2015; O'Rourke et al., 2018) (Table 1).

Table 1

Main parameters for in vitro experiments

	Experimental values
Cell seeding density	0.5×10^6,0.78×10^6,1.5×10^6 cells/ml
Equivalent cell seeding density after stabilization	20×10^6,31×10^6,60×10^6 cells/ml
Oxygen concentration	1%, 3%, 7%, or 21%
Glucose concentration	25 mM
Duration	24 h

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 phase‐contrast water immersion objective. High‐throughput quantification of cell viability from 3D image stacks by adapting a readily available ImageJ protocol was performed.

The spatial distribution of viable cells within the constructs was also estimated. From the 3D image stacks three different zones that correspond to the top, middle, and bottom of the gels were identified and the mean value of live and dead cells per zone was calculated by analyzing seven stacks per zone.

Glucose detection assay

Glucose consumption was quantified by an enzymatic assay (Glucose [HK] Assay Kit, GAHK20, Sigma Aldrich) according to the manufacturer's guidelines. Briefly, after 24 h incubation of the cellular gels, the supernatants were collected for further analysis. The reconstituted reagent was added to each sample and the resulting solution was incubated for minutes at room temperature. During that time, glucose was phosphorylated by ATP, a reaction which was catalyzed by hexokinase. Glucose‐6‐phosphate (G6P) was then catalytically oxidized to 6‐phospho‐gluconate in the presence of oxidized nicotinamide adenine dinucleotide (NAD). Due to this oxidation, an equimolar amount of NAD was reduced to NADH, thereby changing the optical absorbance of the sample. The consequent increase in absorbance was measured at 340 nm and was directly proportional to glucose concentration.

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 diffusion‐reaction equations that describe the interactions between oxygen (), glucose () and VEGF () concentrations and the cell population (), 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.

Figure 1

Cell culture well schematic.

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). where represents the oxygen concentration in the gel, the local cell density in the gel, the diffusion coefficient of oxygen in the gel, the maximum oxygen consumption rate by the cells, and 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 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 represents the glucose concentration in the gel. Here, is the diffusion coefficient of glucose in the gel, the maximum glucose consumption rate by the cells, and 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 conditions, . The functional form of this new term is based on supplemental consumption of glucose in anaerobic conditions happening around the same oxygen level as the weakening of the oxygen metabolic reaction () but being negligible for high oxygen concentration conditions. 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 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 represents the VEGF concentration in the gel, is the diffusion coefficient of VEGF in the gel, the VEGF degradation rate and 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 is the hypoxic threshold for VEGF production and represents a crowding factor for the cells. VEGF is also assumed to diffuse freely in the medium: where represents the VEGF concentration in the medium and the diffusion coefficient of VEGF in the medium.

Figure 4

Functional analysis of dCTX0E03 cells under different oxygen conditions. (a) Glucose consumption by dCTX0E03 cells seeded in collagen and cultured under a range of ambient oxygen concentrations for 24 h. Glucose concentration was quantified using a biochemical assay. Data expressed as means ± SEM (n = 4 independent repeats, three samples per condition). (b) VEGF release from dCTX0E03 cells seeded in collagen and cultured under a range of ambient oxygen concentrations for 24 h. VEGF concentration was measured via ELISA. 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%).

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 controls the oxygen‐related death, controls the glucose‐related death 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 (Figure 1). Conditions are imposed at the boundaries to capture the geometrical constraints and relevant transport characteristics of the setup. At the boundary between the cellular gel and the medium, we assume continuity of concentration and flux for oxygen, glucose, and VEGF. Zero flux boundary conditions are imposed for oxygen at the bottom and the sides of the well, whereas the concentration of oxygen on the medium‐air interface was set as constant and equal to the ambient oxygen as prescribed during the experiments. For VEGF and glucose, zero flux conditions were set at the bottom, sides and air interface. The initial oxygen concentration in the gel was set at due to preconditioning of the cell culture medium. Based on the experimental data, the initial glucose concentration at the medium was set as and in the gel, the initial VEGF concentration was set at , and the initial cell density was varied to match the experimentally‐imposed seeding conditions.

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.

Table 2

Initial parameter range

	Nominal value and approximated bounds
Cell density parameters
Maximal cell density (n_max)	n_max=77×10^6 cell/ml
Proliferation rate constant	0
Baseline cell death rate constant (δ_0g)	[3.3 × 10−7, 1.1334 × 10−5 1/s] (Chung et al., 2006; Coy et al., 2020)
Hypoxic cell death rate constant (δ_c)	a[1.7 × 10−8, 4 × 10−6 ]1/s
Glucose deprivation‐induced death rate constant (δ_g)	a[1.7 × 10−8, 4 × 10−6 ]1/s
Oxygen concentration parameters
Diffusion coefficient for oxygen in medium (D_c_m)	D_c_m= 2.62 × 10−5 [1 × 10−5, 4 × 10−5] cm2/s (Han & Bartels, 1996)
Diffusion coefficient for oxygen in gel (D_c_g)	[1 × 10−6, 4 × 10−6] cm2/s (Cheema et al., 2012)
Concentration at which oxygen consumption is ½ maximal (c®)	[6.66×10^−9,4×10^−8] mol/ml (Coy et al., 2020) (0.5% O2)
Maximal rate of oxygen consumption (M_c)	b[1×10^−18,7.7×10^−16] mol/cell/s (Herculano‐Houzel, 2011; McMurtrey, 2015; Streeter & Cheema, 2011; Wagner et al., 2011)
Glucose concentration parameters
Diffusion coefficient for glucose in medium (D_s_m)	[5.65 ×10^−6, 1.09 ×10^−5] cm2/s (Shipley et al., 2009; Suhaimi et al., 2015)
Diffusion coefficient for glucose in gel (D_s_g)	[0.23 ×10^−6, 1.51 ×10^−6] cm2/s (Cochran et al., 2006; Wu et al., 2005)
Concentration at which glucose consumption is 1/2 maximal (s®)	[2,10] mM (Barros et al., 2007; Duarte et al., 2009; Van Zijl et al., 1997)
Maximal rate of Glucose consumption (M_s)	[5 × 10−17, 2.2 × 10−16] mol/cell/s (Gu et al., 2016; McMurtrey, 2015)
Anaerobic threshold for glucose consumption (A)	N/A
VEGF concentration parameters
Diffusion coefficient for VEGF in medium (D_v_m)	[1.3×10^−6,2×10^−6] cm2/s (Mac Gabhann et al., 2005; Mac Gabhann et al., 2007)
Diffusion coefficient for VEGF in gel (D_vg)	[2.9 × 10^−7, 1.13 × 10^−6] cm2/s (Chen et al., 2007; Köhn‐Luque et al., 2013; Wang et al., 2020)
VEGF degradation rate (K)	[2.67 × 10−6, 1.28 × 10−4] 1/s
Hypoxia threshold for VEGF secretion (c_τ)	a[0.1, 2] %O2
VEGF crowding factor (n_τ)	a[1, 60] ×10^6 cell/ml
Baseline VEGF secretion rate at low oxygen (a)	N/A
Baseline VEGF secretion rate based on oxygen (β)	N/A

Based on experimental observations.

Bounds of parameter were adjusted to account for decreased oxygen consumption rate in 3D culture systems (Magliaro et al., 2019).

2.62 10−5

Sensitivity analysis

The model defined by Equations (1)–(9), while built with minimal components, still includes 19 parameters, excluding initial and boundary conditions. The large number of model parameters compared to the relatively small data set, due to the limited spatial and temporal resolution of the in vitro model, leads to an underdetermined system. To help regularize the problem, we performed a sensitivity analysis to prioritize the importance of parameters in predicting a relevant quantity of interest, chosen here as the average concentrations and cell density in the gel, where represents the local concentration of the different species in the gel and the volume of the gel. By combining this prioritization exercise with existing knowledge of the different parameters, we can define realistic intervals for each parameter value. Such intervals are then explored during the optimization procedure (small intervals for impactful parameters, large intervals for minor parameters) by comparing model predictions against experimental measurements.

Multiple approaches exist to perform sensitivity analysis (Saltelli et al., 2004). Given the size of the parameter set we selected the Morris Screening Method (Morris, 1991) which qualitatively evaluates the global sensitivity of each parameter, including coupling and nonlinearity, using a set of elementary effects. More specifically, we use the implementation proposed by Campolongo et al., who developed an optimized parameter sampling algorithm for sensitivity analysis with decreased computational cost (Campolongo et al., 2007).

Input values of the variables of interest are determined based on a sampling algorithm which starts at randomly selected points in the k‐dimensional space and creates a trajectory through all the k‐dimensions. First, the elementary effects (EE i = 1, …, r, j = 1, …, k) used for the Morris screening test are individually computed for each trajectory and each variable of interest.

An elementary effect can be computed if and only if is still in the parameter range. This means that there are elementary contributions. E values are then used to calculate the final sensitivity measures such as the mean absolute value of the elementary effect and the standard. More specifically, and with 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, the set or subset of species of interest (i.e., oxygen, glucose, VEGF, and cell density) and where 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 is the size of the experimental set for species , the vector containing the corresponding predicted values for a given vector of parameter value and vector of initial and boundary values , and 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 one‐dimensional. As for the axial direction, we devised a nonuniform two‐part 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.

Based on a mesh convergence analysis, mesh cells (80% in the medium and 20% in the gel) with a timestep enabled gradient fields to be sufficiently resolved in both domains. Further increasing the spatial or temporal resolution resulted in at most in ~1% change for the average concentration in the gel after 24 h. Next, the Morris sensitivity analysis was run using the open source library SALib (Herman & Usher, 2017) using the intervals presented in 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) we compute the average oxygen in the gel every 0.5 h for 24 h for: (a) ambient oxygen concentrations 1%, 3%, and 7% and for an initial cell density  60 (, (b) average glucose concentration in the culture medium after 24 h for 1%, 3%, 7%, and 19% and for initial cell densities 20, 31, 60 (, (c) average cell density in the gel after 24 h for 1%, 3%, 7%, and 19% ambient oxygen concentrations and for initial cell densities  20, 31, 60 (. For consistency in calculating the cost function (Equation 11), all results were nondimensionalised using the initial concentration for each species. Similarly, for step (2), we compute the average VEGF concentration in the culture medium for ambient oxygen concentrations 1%, 3%, 7%, and 19% and for initial cell densities 20, 31, 60 (.

RESULTS

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% O2 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.

Figure 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 × 106 cells/ml density after stabilization). Data expressed as means ± SEM (n = 4 independent repeats, three samples per condition).

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).

Oxygen consumption characteristics in 3D constructs

Figure 3 displays the temporal changes in the oxygen concentrations measured at the center of cellular stabilized collagen constructs cultured under different ambient oxygen levels. Studies on acellular constructs (Figure 3a) show that the oxygen concentration in the gel equilibrated to ambient levels within 5 h. Any differences in these profiles for cellular gels must be due to cellular metabolism.

Figure 3

Oxygen levels in the center of (a) acellular (n = 1) or dCTX0E03‐seeded constructs (60 × 106 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).

Cellular constructs exhibited time‐dependent oxygen depletion in their core (Figure 3b–d). There was a rapid fall of oxygen toward approximately steady‐state values, with the rates being affected by the ambient oxygen level. The lowest ambient oxygen concentration of 1% caused the steeper gradients for collagen gels with dCTX0E03 cells. Interestingly, the oxygen concentration appears to reincrease after 12 h for cellular constructs cultured at 7%. Between 0 and 12 h, 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 and  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).

Figure 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.

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 oxygen‐related death rate, and the baseline death rate (Figure 5a). The maximal rate of oxygen consumption has the strongest influence on model predictions and interacts most with other parameters. For the VEGF concentration, most parameters exhibit a nonlinear influence on and/or interactions with other parameters (σ/μ* > 0.5). The final concentration is mostly affected by the crowding factor and VEGF degradation (Figure 5b). The effect of a large group of parameters including VEGF production rate, maximum oxygen consumption rate, and baseline death rate can be considered secondary but nonnegligible. For glucose (Figure 5c), the model output varies nonlinearly but monotonically or almost monotonically with the maximum rate of glucose consumption, anaerobic threshold, glucose diffusion in the medium, and the glucose concentration for which consumption is half of the maximum. The global sensitivity analysis also indicates that the maximum rate of oxygen consumption is a potentially influential parameter in glucose concentration. Finally, the cell viability after 24 h is mostly affected by the baseline death rate (Figure 5d). Those parameters identified as having minimal influence on model outputs were then fixed at the nominal values provided in Table 2 ().

Next, we used a particle swarm algorithm to minimize the difference between model predictions and experimental data, via the choice of the remaining parameters , as described in the Methods. Table 3 summarizes the final set of parameters found.

Table 3

Final parameter values

	Final value (practical units)	Final value (modeling units)
Cell density parameters
Maximal cell density (n_max)	n_max=6.0×10^7 cell/ml	n_max=6.0×10^13 cell/m^3
Proliferation rate constant	0	0
Baseline cell death rate constant (δ_0g)	δ_0g=3.18×10^−6 1/s	δ_0g=3.18×10^−6 1/s
Hypoxic cell death rate constant (δ_c)	δ_c=2.53×10^−6 1/s	δ_c=2.53×10^−6 1/s
Glucose deprivation induced death rate constant (δ_g)	δ_g=5.6×10^−7 1/s	δ_g=5.6×10^−7 1/s
Oxygen concentration parameters
Diffusion coefficient for oxygen in medium (D_c_m)	D_c_m=2×10^−9 m2/s	D_c_m=2×10^−9 m2/s
Diffusion coefficient for oxygen in gel (D_c_g)	D_c_g=4.98×10^−10 m2/s	D_c_g=4.98×10^−10 m2/s
Concentration at which oxygen consumption is 1/2 maximal (c®)	c®=1.65×10^−8 mol/ml (1.24% O2)	c¯=5.13 × 10^‐4 kg/m^3
Maximal rate of oxygen consumption (M_c)	M_c=1.84×10^−18 mol/cell/s	M_c=5.88 × 10^‐20  kg/cell/s
Glucose concentration parameters
Diffusion coefficient for glucose in medium (D_s_m)	D_s_m=9×10^−11 m^2/s	D_s_m=9×10^−11  m^2/s
Diffusion coefficient for glucose in gel (D_s_g)	D_s_g=2.67×10^−10 m^2/s	D_s_g=2.67×10^−10 m^2/s
Concentration at which glucose consumption is 1/2 maximal (s®)	s¯=7.8 mM	s¯=1.39 kg/m^3
Maximal rate of glucose consumption (M_s)	M_s=9.8×10^−18 mol/cell/s	M_s=1.75×10^−18 kg/cell/s
Anaerobic threshold for glucose consumption (A)	A=4.6	A=4.6
VEGF concentration parameters
Diffusion coefficient for VEGF in medium (D_v_m)	D_v_m=1.32×10^−10 m2/s	D_v_m=1.32×10^‐10 m2/s
Diffusion coefficient for VEGF in gel (D_vg)	D_vg=4.16×10^−11 m2/s	D_vg=4.16×10^−11 m2/s
VEGF degradation rate (K)	K =8.37×10^−5 1/s	K =8.37×10^−5 1/s
Hypoxia threshold for VEGF secretion (c_τ)	c_τ=1.40×10^−8 mol/ml (1.08% O2)	c_τ=4.49×10^‐4   kg/m^3
VEGF crowding factor (n_τ)	n_τ=2.32×10^7 cell/ml	n_τ=2.32×10^13   cell/m^‐3
Baseline VEGF secretion rate at low oxygen (a)	α=9.24×10^−2 pg/cell/s (molO2/ml)−1 (1.21×10^−9 pg/cell/s/% O2)	α=2.92×10^−21   kg/cell/s (kgO2/m3)−1
Baseline VEGF secretion rate based on oxygen (β)	β=2.86×10^−7 pg/cell/s	β=2.86×10^−22   kg/cell/s

Figures 6 and 7 compare simulation predictions and experimental data for the final optimized parameter value set. Overall, the model replicates the general trends for the viable cell density, nutrient consumption, and VEGF release. The best fit for cell viability is for the 7% O2 data set. For the remaining conditions, 1% and 3% O2, the model tends to respectively underestimate and overestimate the mean viable cell density. With regard to glucose consumption, the fit against experimental observations appears to be worse for  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  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.

Figure 6

Comparison between experimental and simulation results based on final parameter values of (a–d) cell viability, (e–h) glucose, and (i–l) VEGF in the media after 24 h.

Figure 7

Comparison between experimental and simulation results based on oxygen profiles in the gel after 24 h.

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 fact that cells adapt by recalibrating their metabolic profile and activating antiapoptotic pathways. Indeed, we observed an increase in the rate of glucose utilization under low O2 tension.

Oxygen and glucose deprivation have also been correlated with changes in growth factor release (Mac Gabhann et al., 2007). Our results are consistent with reports that VEGF expression increases under hypoxic conditions; a response that has been linked to neuronal protection and nerve regeneration (Jin et al., 2001; 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 computationally effective model to describe interactions between dCTX0E03 cells and soluble factors that can be readily used to investigate various nerve repair scenarios. The functional forms of the equations were adapted from previous cell‐solute mathematical models (Chung et al., 2006; Coy et al., 2020; McMurtrey, 2015; Streeter & Cheema, 2011) developed for other cell types. Much of the modeling work done in tissue engineering so far has produced interesting results and generated general hypotheses about the optimization of tissue‐engineered constructs or the tissue culture conditions. However, many of the mathematical models were not benchmarked against a specific or consistent experimental data set or were validated by comparing theoretical simulations with scattered data from multiple sources from the literature, often for different cell types. Here, our focus was to parametrize the model using dedicated in vitro experiments. Aside from the novel set of differential equations that make up the model, we also optimized the parameter values and tailored them to the metabolic and functional characteristics of dCTX0E03s within collagen constructs. The final derived values are mostly within the range reported in the literature (Table 2) for other cell types. One noteworthy exception is the maximal rate of glucose consumption , 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  cells/ml at 3%, 7% O2. 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 this study, only the effect of oxygen was included in the VEGF governing equation.

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, 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.

CONFLICTS OF INTEREST

The authors declare no conflicts of interest.