Skip to main content
Advertisement
  • Loading metrics

Modelling-informed cell-seeded nerve repair construct designs for treating peripheral nerve injuries

Abstract

Millions of people worldwide are affected by peripheral nerve injuries (PNI), involving billions of dollars in healthcare costs. Common outcomes for patients include paralysis and loss of sensation, often leading to lifelong pain and disability. Engineered Neural Tissue (EngNT) is being developed as an alternative to the current treatments for large-gap PNIs that show underwhelming functional recovery in many cases. EngNT repair constructs are composed of a stabilised hydrogel cylinder, surrounded by a sheath of material, to mimic the properties of nerve tissue. The technology also enables the spatial seeding of therapeutic cells in the hydrogel to promote nerve regeneration. The identification of mechanisms leading to maximal nerve regeneration and to functional recovery is a central challenge in the design of EngNT repair constructs. Using in vivo experiments in isolation is costly and time-consuming, offering a limited insight on the mechanisms underlying the performance of a given repair construct. To bridge this gap, we derive a cell-solute model and apply it to the case of EngNT repair constructs seeded with therapeutic cells which produce vascular endothelial growth factor (VEGF) under low oxygen conditions to promote vascularisation in the construct. The model comprises a set of coupled non-linear diffusion-reaction equations describing the evolving cell population along with its interactions with oxygen and VEGF fields during the first 24h after transplant into the nerve injury site. This model allows us to evaluate a wide range of repair construct designs (e.g. cell-seeding strategy, sheath material, culture conditions), the idea being that designs performing well over a short timescale could be shortlisted for in vivo trials. In particular, our results suggest that seeding cells beyond a certain density threshold is detrimental regardless of the situation considered, opening new avenues for future nerve tissue engineering.

Author summary

Mathematical models are increasingly gaining traction in biomedical sciences in general and in peripheral nerve repair in particular. These models can both help to unveil mechanisms underlying nerve regeneration post-injury and inform the design of new repair strategies, while continuously being improved through experiments and clinical data. In particular, cell-solute models provide flexible, computationally-efficient frameworks that enable access to spatio-temporal information difficult to reach in vivo in nerve repair scenarios. In this work, we derive such a model and apply it to the case of Engineered Neural Tissue repair constructs after their transplant into the nerve injury site. One key feature of such repair constructs is the spatial seeding of therapeutic cells capable of producing signalling molecules that promote nerve regeneration. We simulate the activity of such cells and evaluate, after 24h, cell survival and signalling molecule output for a range of nerve repair construct designs. The underlying hypothesis is that designs performing well at such short timescales could be shortlisted for in vivo trials, the latter being costly and time-consuming. Our results suggest that seeding cells beyond a certain density threshold is detrimental regardless of the situation considered, opening new avenues for future nerve repair construct design.

Introduction

More than one million people per year are affected by peripheral nerve injuries (PNI) in Europe and the USA [1] and healthcare costs are estimated to exceed £1Bn per year in the USA alone [2]. Paralysis and loss of sensation are hallmarks of severe PNIs, which can lead to lifelong pain and disability for patients. Neurons in the peripheral nervous system have some ability to regenerate following damage, via a process involving both neuronal and non-neuronal cells and signalling molecules.

The extent of nerve regeneration is largely determined by the severity of the injury, which often requires a surgical intervention. Gold-standard treatments for large-gap PNIs involve an autograft, i.e. the transplantation of healthy nerve sections from elsewhere in a patient’s body into the injury-induced nerve gap. However, studies have reported that only around 50% of autograft patients regain a good level of function [36]. In addition, nerve and tissue damage at the donor nerve site results in scarring and loss of function and can lead to more pain and increased risk of infection. The need to harvest a donor nerve also lengthens the time spent in surgery and the number of eligible nerve donor sites is limited.

Bioengineered peripheral Nerve Repair Constructs (NRCs) are being developed as alternatives to address these issues. A large range of synthetic and natural biomaterial bases for these constructs have been explored [7,8]; here we focus on one example with high translational potential, Engineered Neural Tissue (EngNT) [9]. EngNT constructs are composed of stabilised hydrogel, typically collagen, containing aligned therapeutic cells. Many cell types have been explored dependent on the therapeutic scenario [10] including differentiated adipose-derived stem cells (dADSC) that are able to produce Vascular Endothelial Growth factor (VEGF) under low oxygen conditions [11,12] to promote revascularisation of the injury site. EngNT constructs mimic the geometry of the host nerve, so are typically cylindrical in shape, and surrounded by a protective sheath of material that can be made permeable to exchange of molecules. These NRCs provide a supportive microenvironment to promote neurite and blood vessel regrowth between nerve stumps, and the technology potentially enables careful spatial seeding of therapeutic cells to facilitate this [13].

In this context, the goal is then to determine the design of NRCs that would yield maximal nerve regeneration and lead to functional recovery. So far, most of the investigations on this question have generally focused on experimentally measuring the benefit of a given NRC design after several days to several weeks after transplant (e.g. therapeutic cell type, sheath material, see references in [10,14]), which is costly and time-consuming. On the other hand, short-term impacts of cell-seeding strategies, culture conditions (e.g. initial concentration of oxygen) or sheath material (e.g. sheath porosity) on cell survival and signalling factor production, which are likely to affect long-term regeneration, are challenging to assess experimentally in vivo because there are so many possible options to test.

Mathematical models have been widely used in the field of regenerative medicine [15], including peripheral nerve repair [16]. They generally provide a description of the interactions between cells and their environment at different scales, including the surrounding substrate, nutrients and signalling factors. For problems that require a description of subcellular mechanisms, discrete models such as Agent-Based Model [17] have been developed to retrieve collective cell behaviours, including tumour [18] and soft tissue growth [19].

For problems requiring a description at the scale of the whole tissue, models based on multiphase theory in porous media [20,21] have been developed. These models consider the cells and their surrounding microenvironment to be a continuum with effective properties. In particular, these models have been used to describe situations analogous to EngNT constructs in the context of tissue culture bioreactors [2227].

Broadly speaking, these models comprise a set of effective, non-linear coupled partial differential equations describing both the spatial and the temporal behaviours of cell populations, solute transport and the porous media formed by the substrate (e.g., structured scaffold, hydrogel). For instance, Lemon et al.[28] developed a 3-phase model (i.e., cell, culture media, porous scaffold) with a focus on mechanical interactions, allowing them to describe cell migration in an engineered tissue. This model was later refined to include cell proliferation and nutrient transport [29,30] making it a general framework to study engineered tissue. Similarly, Ochoa et al.[31] used formal Volume Averaging [21], to derive a model of solute transport in an idealized cellular media and were able to obtain representative transport equations at the scale of the whole media. Such a model was then extended to biofilms [32] and to culture bioreactors [33,34], taking into account cell metabolism and cell proliferation, making it a strong candidate for cell-solute descriptions for engineered tissue.

Continuous models can therefore be seen as valuable tools to study engineered tissue at a limited cost. That is why in the present work, given the typical size of NRC constructs (>1cm), the number of cells usually seeded (between 3×104 and 1.2×106 [35,36]) and the time scale we are interested in (24h), we derive a continuous model focusing on cell-solute interactions in an EngNT seeded with dADSCs, since the latter are relevant for NRCs currently under development and have been characterised in vitro [12]. We point out that, even though we focus on specific types of constructs and therapeutic cells, we derive the cell-solute model in a relatively general manner that can easily be extended to other repair scenarios. The objective of such a model is then to assess the performance of different NRC designs within the first 24h following transplant.

Briefly, during this period of time the therapeutic cells proliferate in the collagen, consuming the oxygen present in their environment. Then, when the oxygen concentration reaches low values, the cells become stressed and produce VEGF that will later promote the migration of endothelial cells and the creation of a vascular network that will perfuse the construct with oxygen and nutrients, supporting both the therapeutic cell population and the subsequent neurite regrowth. During the first 24h, we therefore focus on simulating the evolving cell population, oxygen, and VEGF concentrations, along with their interactions, before endothelial cell migration occurs. We then determine which design yields the maximal mean cell density, as therapeutic cells, besides VEGF secretion, can also provide mechanical support via the production of extra cellular matrix proteins [37,38] and directional cues [9] for vessel sprouts and neurites. We also determine which design yields the maximal mean VEGF concentration and maximal VEGF gradients after 24h, as VEGF plays a central role in regenerative angiogenesis, which is involved in subsequent stages of nerve repair [39].

In doing so, the underlying ideas are that 1) designs that perform poorly on this short-time scale could be discarded, hence saving experimental time and 2) modelling could help us understand the underlying mechanisms influencing performance and could indicate new avenues for design improvement.

The model itself is parametrised using in vitro data when possible and left uncertain otherwise. Doing so allows us to evaluate in a systematic manner a wide range of NRC designs. In particular, we investigate the impact of different cell-seeding strategies (initial cell density and initial spatial distribution of cells), as therapeutic cell distribution is likely to control VEGF secretion and cell survival. We also investigate the influence of sheath material (e.g., porosity, thickness), culture conditions pre-implantation (e.g. initial oxygen concentration) and surrounding tissue conditions, on cell survival and vascular growth factor production.

Methods

Nerve repair construct geometry

We consider the NRC to be composed of a cell-seeded collagen cylindrical core of length L (m) and radius R1 (m) surrounded by an annular, acellular sheath of materials of thickness T (m), as illustrated in Fig 1. This geometry mimics a typical repair scenario in an in vivo test model, such as the rat sciatic nerve. The sheath itself can be impermeable (Fig 1 when Lp = 0), partially or entirely (Fig 1 when 2Lp = L) permeable to molecular exchanges with the surrounding tissue. We further assume that the core and the sheath form two distinct continuous media with properties that can be described using an effective description.

thumbnail
Fig 1. Nerve repair construct geometry.

Sheath partially made of porous material permeable to molecular exchanges at the construct extremities. The parameter R1 denotes the radius of the cell-seeded collagen cylindrical core, T the thickness of the sheath, R2 = R1+T the radius of the nerve repair construct, L the length of the construct and Lp the length of the porous section of the sheath (2Lp in total).

https://doi.org/10.1371/journal.pcbi.1009142.g001

Cell density dynamics

Different therapeutic cell types have been used in EngNT designs [13,40]. As mentioned in the Introduction, we focus on the case of dADSCs because they are relevant for EngNT NRCs currently under development and have been characterised in vitro. However, the cell-solute model is derived in a relatively generic manner and could be readily adapted for other cell types.

In general, the local cell density is determined by the growth and decay rates of the cell population, as well as cell migration. In stabilised collagen gels, however, cell migration can take weeks, and in vitro experiments using different cells in similar materials have shown no significant cell displacements over the first five days following the seeding [41]. Given that we focus on the first 24 hours post implantation, in this work we neglect cell migration mechanisms and write the following balance for the local cell density (1) where n represents the local cell density (cell∙m-3), Gn the cell density growth (cell∙m-3∙s-1), and Qn the cell density decay (cell∙m-3∙s-1). Following the work of Coy et al. [12], we assume that cell density growth has the following form (2) where cg denotes the local oxygen concentration in the gel (kg∙m-3), β the cell growth rate (m3∙kg-1∙s-1) and nmax the maximum cell density (cell∙m-3). We note that the cell growth term depends only on local cell density and oxygen concentration since we have assumed all other substrates, such as glucose, to be in excess. We also note that the local cell density follows a logistic growth form to account for nutrient and space competition, hence the inclusion of nmax. With regards to cell decay we also follow the work of Coy et al. [12] who assumed a linear decay term, yielding (3) where δ denotes the cell death rate (s-1). We note that δ is in fact an increasing linear function of the initial cell density to represent the effects of initial crowding on later time points and provide its expression in the parameter value section below.

Molecular transport in the cell-seeded collagen gel

PNI regeneration involves a large range of molecules that interact with each other and the embedded cells. In this work, however, we focus on oxygen and vascular growth factors, as these molecules along with their interactions are central in cell survival and angiogenesis. Other nutrients such as glucose are assumed to be in excess for the first 24h, and the impact of other growth factors such as brain derived neurotrophic factor, glial cell line-derived neurotrophic factor or nerve growth factor are neglected. However, the framework we establish is sufficiently general that these molecules could be included in future studies.

In this context, we consider the transport of oxygen to be driven by molecular diffusion and assume that cells metabolize oxygen following Michaelis-Menten kinetics, as is commonly used in the literature for oxygen-limited conditions [4245], so that (4) where cg corresponds to the local oxygen concentration (kg∙m-3), Dc,g to the diffusion coefficient of oxygen within the collagen gel (m2∙s-1), M to the maximal metabolic rate (kg∙cell-1∙s-1) and c1/2 to the oxygen concentration (kg∙m-3) for which the metabolic reaction equals half its maximal value.

Vascular growth factors form a family of molecules. In this work, however, we focus on modelling vascular endothelial growth factor-A (VEGF-A), which is regarded as the most influential molecule of the VEGF family upon angiogenesis [46,47]. In the following, VEGF-A is simply referred to as VEGF for simplicity.

In this context we assume the transport of VEGF in the cell-seeded collagen gel to be driven by molecular diffusion. Contrarily to oxygen, which is a stable molecule, VEGF, which is secreted by the cells, is an unstable protein and its degradation is usually modelled using a linear decay term [4852] so that (5) where vg corresponds to the local VEGF concentration (kg∙m-3), Dv,g to the diffusion coefficient of VEGF within the collagen gel (m2∙s-1), K to the VEGF degradation rate (s−1) and Gv(n,c) to the production rate (kg∙m-3∙s-1) of VEGF. With regards to VEGF production, we use the expression derived by Coy et al. [12], which takes into account upregulation of VEGF in low-oxygen transport regimes so that (6) where α denotes the baseline VEGF production rate (kg∙cell-1∙s-1), ch the oxygen threshold under which VEGF production is upregulated by oxygen (kg∙m-3), Vm the factor controlling the magnitude of such a regulation and B (m3∙kg-1) the factor controlling the transition between baseline and regulated states. We point out that both α and Vm are in fact functions of the initial cell density, and that their expressions are given in the parameter value section below.

Finally, we note that, in general, the diffusion coefficients in Eqs 4 and 5 should depend on the cell density, as an increase in cell density implies a higher tortuosity in the gel and anomalous diffusion due to interfacial effects. Practically however, this yields a significant impact only for high cell densities [44,53], which are not the focus of this work. Analogous to cell migration mechanisms, we also neglect effects of collagen matrix remodelling on molecular transport, assuming that it takes place on a longer timescale. Hence, we assume the diffusion coefficients to be independent of the cell density in the gel and the evolving collagen structure.

Molecular transport in the annular sheath

We consider the sheath surrounding the gel to be a porous media devoid of any cells and its porosity to be a free parameter ranging from 0 (completely impermeable to molecular exchanges) to 1 (completely permeable). We further consider the pores to be larger than the characteristic size of the molecules transported. In this situation we can no longer neglect a priori the effects of porosity and tortuosity on molecular transport [54]. Therefore, the transport equations for oxygen and VEGF in the sheath are given by (7) for oxygen and (8) for VEGF, where Dc,s and Dv,s denote the molecular diffusion coefficient of oxygen and VEGF respectively (m2∙s-1) in the sheath, ϵ denotes the sheath porosity and τ its tortuosity. We recall that the tortuosity is a function of the underlying porous structure, which varies widely according to the materials and methods used to create the NRC sheath [5557]. In this work we propose to take a step back and narrow our focus to archetypal structures, commonly encountered in porous media, for which it is possible to derive theoretically an expression for the tortuosity. Precisely, we will consider the cases of overlapping spheres [58], cylinders [59] and disordered tubes [60]. Table 1 summarizes the different expressions for the tortuosity depending on the aforementioned porous structures.

thumbnail
Table 1. Tortuosity expressions as a function of porous structures.

https://doi.org/10.1371/journal.pcbi.1009142.t001

Initial and boundary conditions

Starting with the envelope of the NRC, we expect the tissue surrounding it to maintain, at least during the first days after the graft, a form of homeostasis and therefore assume that both oxygen and VEGF fields are constant outside of the construct. We note that, practically, concentration values at the severed nerve ends and in the surrounding tissue are likely to be different. Due to the lack of measures however, we set them at a value, representative of nerve tissue in general. This translates into the following Dirichlet conditions (9) (10) for the oxygen and VEGF concentration fields associated to the gel and (11) (12) (13) (14) for the oxygen and VEGF concentration fields associated to the annular sheath.

Next, with regards to the internal boundaries, we impose conservation of the molecular flux at the interface between the cylinder core and the sheath (r = R1) for both the oxygen and VEGF fields (15) (16) In addition, we assume that the diffusion process is slow enough so that the local thermodynamic equilibrium is reached at the interface between gel and sheath material yielding (17) (18) where λc and λv denote the partition coefficients between the gel and the sheath for the oxygen and VEGF respectively. We note that there is no need for boundary conditions for the cell density since we assumed no migration mechanisms.

With regards to initial conditions, we prescribe a uniform concentration throughout the construct for both oxygen and VEGF to represent storage conditions before implantation into the injury site, so that (19) (20) Finally, we impose the initial cell density to be uniform across the cylinder core cross-section but still able to depend upon axial position, yielding (21) leaving the opportunity to explore the consequences of non-uniform initial cell-seeding on oxygen and VEGF distributions.

Cell seeding strategy evaluation

One of the main aims of this work is to identify cell seeding strategies that could increase cell density, along with VEGF concentrations and gradients within NRCs that may favour vascular regrowth. To that end, we first define the following quantities in the cylindrical core (22) (23) (24) (25) where denotes the average cell density, the average oxygen concentration, the average VEGF concentration over V, the volume of the collagen core. vSD represents the standard deviation associated to the VEGF concentration through the collagen core. vSD is used here an indirect measure of VEGF gradients. We then define initial cell densities that will maximise and vSD the following way (26) (27) (28) where and represent the initial cell densities leading to the maximal average cell density, maximal average VEGF concentration and maximal VEGF spreading at time t respectively.

Simulations

Simulations of different cell densities and distributions over a NRC geometry along with prediction of oxygen and VEGF concentrations are carried out using finite element methods using COMSOL Multiphysics.

Parameter values

The parameters associated with this problem can be divided into two classes. First there are the parameters that have been previously determined using experiments, such as those associated with the transport of oxygen and VEGF in the collagen gel [12]. Then, there are parameters for which values are rather taken from a range, whether due to lack of contextual information (e.g. VEGF and oxygen concentrations in the tissue surrounding the NRC) or due to deliberate design choices (e.g. initial cell densities, sheath thickness). In the following tables, parameters are given defined values when possible and ranges of values otherwise. We also point out that Tables 26 present parameter values in two different systems of units; practical units, i.e. units useful for discussing and disseminating results (e.g. concentrations in %O2, ng∙ml-1), and modelling units, i.e. SI units useful for solving mass balance equations (e.g. concentrations in kg∙m-3).

Starting with the NRC itself, Table 2 summarizes the characteristics of the cylinder core and annular sheath geometry. We note that the choice for the NRC length L and cylinder radius R1 was made according to NRC designs currently under development. The range of values for the sheath thickness T was made to reflect existing designs [55,6166].

Table 3 presents parameters associated with the cell population dynamics (Eqs 13). We recall that the cell death rate δ is a function of the initial cell density and has the following expression [12] (29) where δ0 and δ1 are fitting parameters which values are reported in Table 3.

Tables 4 and 5 summarize the parameter values associated with the transport of oxygen and VEGF in the gel and the sheath (Eqs 48). We note that we consider the phase filling the pores of the sheath to be similar to water, and hence choose the value for the diffusion coefficient accordingly. We also recall that the baseline VEGF production rate α and the regulation magnitude factor Vm are functions of the initial cell density n0 and have the following expressions [12] (30) (31) where α0, α1, α2, Vm,0 and Vm,1 are fitting parameters with values are reported in Table 5.

Finally, Table 6 summarizes the initial and boundary condition values. We point out that due to the lack of oxygen concentration data available for the rat sciatic nerve model, the choice for the range of tissue oxygen concentration ctissue was made to encompass the physiological oxygen range of 4–7% O2 measured in vivo in humans [70]. Similarly, the range of values for tissue VEGF concentration is based on measurements in serum (1.70×10-2 to 3.00×10-1 ng∙ml-1) and plasma (1.70×10-2 to 3.00*10-1 ng∙ml-1) of healthy controls [71]. We also note that the initial oxygen concentration ranges from 1% to 21%O2 to reflect possible storage conditions pre-implantation into the injury site, and that initial VEGF concentration ranges similarly to the tissue VEGF concentration for consistency. As for initial cell density, the range of values used encompasses the average in vivo Schwan cell density in normal nerve (2.00×107 cell∙ml-1)[35,72], as well as typical cell densities used in NRCs (2.00×107 to 1.60×108 cell∙ml-1) [35,36]. Finally, following the work of Coy et al. [12] we assume that the partition coefficients λc and λv are close to one so that Eqs 17 and 18 simplify to continuity of the concentrations between gel and sheath.

Results

We simulate cell-solute interactions in an array of nerve repair scenarios using the model derived in the method section. Each scenario has a different setting but the same structure; we seed an increasing number of cells in the collagen gel and determine which initial cell density yields the highest mean cell density (), highest mean VEGF concentration () and largest standard deviation in VEGF concentration () 24h after transplant. We recall that is used as an indirect measure of VEGF gradient at the scale of the NRC. As mentioned in the Introduction, mean cell density, VEGF concentration and VEGF gradients are involved in regenerative angiogenesis which is one of the subsequent steps in nerve repair.

The scenarios are organised as follows. We start by assuming the sheath surrounding the collagen core to be impermeable and explore a range of storage conditions (i.e. initial oxygen and VEGF concentrations) and tissue conditions (i.e. oxygen and VEGF concentrations at the nerve stumps). We then estimate the effects of sheath porosity, sheath material (i.e. the underlying porous structure) and sheath thickness. Finally, we investigate the cases of partially porous sheaths (i.e. sheath with different porous length) and non-uniform seeding strategies (e.g. seeding cell preferentially in the centre of the construct or close to the nerve stumps).

In each scenario, unless specified, the initial concentration of oxygen in the collagen gel and sheath is equal to 21%O2, the surrounding tissue concentration of oxygen is equal to 5%O2 and both the initial and tissue VEGF concentrations are equal to zero. Further, when applicable, the sheath is assumed to have a globular structure that could be approximated by a body of overlapping spheres with an associated porosity ϵ = 0.8 and a thickness T = 0.25mm.

Impermeable sheath

In this section we assume the sheath surrounding the gel to be impermeable so that molecular exchanges between the NRC and surrounding tissue only happen at the proximal and distal stumps. The underlying idea here is that an impermeable sheath could prevent the loss of VEGF to the tissue, hence increasing its concentration in the construct, with potential benefits in stimulating blood vessel growth.

Fig 2 shows heatmaps of cell density, VEGF concentration and oxygen concentration over 24h for an initial uniform cell density n0 = 1.78×108 cells/ml. The cell density globally decreases non-uniformly over 24 hours, generally with regions of higher cell density close to proximal and distal stumps. This global decrease is mainly due to the longitudinal diffusion time of oxygen

thumbnail
Fig 2. Cell-solute interactions in an impermeable nerve repair construct result in non-uniform distributions of cells, oxygen and VEGF.

Heatmaps of cell density, oxygen and VEGF concentration over 24 hours for an initial cell density n0 = 1.78×108 cells/ml. Only the collagen gel cylindrical core is displayed as the surrounding sheath is impermeable.

https://doi.org/10.1371/journal.pcbi.1009142.g002

(32) being much longer (tc,L≈34.5h≈1.44 days) than the consumption time of oxygen (tc,M≈1.95h) (33) This means that oxygen is consumed faster than it is delivered from the surrounding tissue. We note that at (Eq 1), so that the cells are actually proliferating. This is due to the high initial concentration of oxygen, which can sustain cell proliferation on a short timescale. After a few hours, however, most of the initial oxygen has been consumed and new oxygen molecules have yet to diffuse from the tissue. This causes a significant decrease in the cell proliferation rate Gn (Eq 2), resulting in the global decrease observed in Fig 2.

With regards to the spatial non-uniformities, for the first 6 hours the cell density is lowest close to proximal and distal stumps. The oxygen concentration in these regions is lower than in the middle of the construct at early times due to the loss of oxygen to the tissue by diffusion (c0>ctissue). This lower concentration comparatively hinders cell proliferation, resulting in an earlier decrease in cell density. After 6h, however, cell densities are the lowest close to the centre of the construct. At this time, cells have consumed most of the oxygen initially present in the construct. Cells close to the proximal and distal stumps, however, have still access to oxygen from the tissue. In addition, cells close to the centre receive only a fraction of the oxygen coming from the tissue as it is first consumed by the cells close to the stumps. These inequalities in oxygen access result in the observed larger cell density close to the extremities and lower cell density in the centre.

Fig 2 also shows that the VEGF concentration increases for the first 12 hours and then decreases in a non-uniform manner. The initial increase is due to the longitudinal diffusion time of VEGF (34) being much longer (tv,L≈139h≈5.8 days) than the secretion time of VEGF (tv,P≈3.63h) defined as (35) where vp is a VEGF concentration (here vp = 100 ng/ml, which corresponds approximately to half the maximum VEGF concentration in the first 24h). Consequently, the VEGF secreted by the cells accumulates inside the construct, hence the increase of VEGF concentration. The long-term decrease is due to the cell density decrease, so that cell secretion of VEGF, despite the oxygen concentration falling below the hypoxic threshold, is no longer enough to compensate for the combined effect of longitudinal diffusion and VEGF degradation.

The VEGF concentration is also spatially non-unform with lower values close to the proximal and distal stumps. This is partly because the concentration of oxygen in the tissue is slightly higher than the hypoxic threshold, causing the cells to secrete less VEGF and partly because of loss to the tissue due to longitudinal diffusion.

Finally, Fig 2 shows that oxygen concentration also decreases non-uniformly. This is mainly due to cell activity, since the time for oxygen metabolism is much smaller (tc,M≈1.95h) axial diffusion time (tc,L≈34.5h≈1.44 days) leading to lower concentration close to the centre of the construct. As for the lower oxygen concentration close to proximal and distal stumps at early times, they are due to loss by diffusion (c0>ctissue).

Whereas Fig 2 focuses on a single initial cell density and does not quantitatively inform on the best cell seeding strategy, Fig 3 explores the impact of initial cell density on cell survival and VEGF concentration after 12h, 24h and 120h (5 days). Specifically, Fig 3A shows the mean (black line) and the range (light grey area) associated with cell density as a function of the initial cell density. The black dots represent the initial cell density value which yielded the largest mean cell density at a given time point, i.e., .

thumbnail
Fig 3. Increasing cell density beyond a threshold decreases both cell survival and VEGF secretion.

A) cell density and B) VEGF concentration means (black lines), standard deviations (dark grey areas) and ranges (light grey areas) as functions of the initial cell density, after 12 hours, 1 day and 5 days. The black dots indicate the initial cell density yielding A) the highest mean cell density, i.e. and B) the highest mean VEGF concentration, i.e. , and the highest standard deviation of VEGF concentration, i.e. .

https://doi.org/10.1371/journal.pcbi.1009142.g003

After 12h, we see that cells/ml, which corresponds to approximately half the maximal initial cell density tested. Consequently, the model predicts that increasing the initial cell density beyond this value deteriorates cell survival. A larger initial cell density implies a faster consumption of the oxygen initially present in the construct (e.g., n0 = 4.00×108 cells/ml yields tc,M = 0.96h). As a consequence, the rate of change of the cell density becomes negative earlier, causing the cell density to decrease sooner. In addition, the cell death rate δ is increasing with n0, causing the high cell density to decrease faster. On the other hand, a lower initial cell density implies a lower cell death rate along with a slower consumption of the oxygen (e.g., n0 = 1.00×107 cells/ml yields tc,M = 38h) and therefore a longer cell proliferation time. We see, however, that this does not necessarily compensate for the impact of the initial gap in cell density, even though seeding n0≤4.00×107cells/ml leads to an increase in the mean cell density after 12h. This balance between initial cell density, cell proliferation and death rate explain the maximum value in the mean cell density after 12h.

Fig 3A also shows that, after 12h, the mean cell density is generally much closer to the lower bound of the cell density range, except for very low initial cell density. This is a direct effect of the non-uniformities discussed in Fig 2. For high initial cell densities, the initial oxygen is consumed quickly and oxygen from the tissue cannot penetrate far enough in the construct due to slow axial diffusion and cell consumption close to the stumps. It is possible to estimate the characteristic penetration length of oxygen (36) where nmax refers to the upper bound of the cell density range which represents, for high initial cell density, the cell density close to the stumps. Such a length is approximately equal to 1.3mm for nmax = 8.50×107cell/ml (i.e. the upper bound associated to after 12h), so that the region perfused with oxygen (z<lc;z>Llc) is significantly smaller than the region deprived of oxygen (lc<z<Llc), explaining why the mean cell density is closer to the lower bound of the cell density range.

By comparison, for a lower initial cell density, the consumption of oxygen is slower, so that the cell density remains higher close to the centre, resulting in the mean cell density being closer to the upper bound of the range.

Finally, Fig 3A shows that after 12h, the cell density range is minimal for n0≈4.80×107 cells/ml and maximal for n0≈2.15×108 cells/ml. This minimal size implies almost uniform cell density. This happens when the region close to the centre starts exhibiting lower cell density than the extremities due to oxygen deprivation. The maximal cell density range, similar to the maximum in mean cell density, is the result of the balance between effects of initial cell density, cell proliferation and cell death rate.

After 24h, Fig 3A shows that the mean cell density has generally decreased but still exhibits the same variations, with a maximum for cells/ml. This decrease is more pronounced for higher initial cell densities due to the continuing lack of oxygen at the centre of the construct and the cell death rate δ being an increasing function of n0. By comparison, lower initial cell densities consume oxygen slower and exhibit lower cell death rates, leading to a smaller decrease over 24h, closing the gap in initial cell density, explaining the decrease of .

Fig 3A also shows that after 24h the boundaries of the range have globally decreased except for the lower bound, which has increased for very low initial cell densities (n0≤2.0×107 cells/ml). This increase being located before the range reaches its minimum size means that cells close to the stumps are still actually proliferating after 24h. This is because the cells in these regions still have access to oxygen from the tissue and are at sufficiently low density so that they don’t have to compete for it.

After 5 days, Fig 3A shows a major change in the cell density behaviour. The cell density reaches zero for n0>2.00×108 cells/ml, with no cells left alive, even close to the stumps. This is a consequence of the parameter values used to solve Eqs 13. To show this, we first determine the stationary states associated with the cell density, which are defined as (37) where nstat and cg,stat denote the stationary cell density and oxygen concentration respectively. Using the expressions for Gn and Qn (Eqs 2 and 3) the above equation has two solutions (38) (39) Performing a stability analysis of the above two solutions [73] yields (40) (41) which means that the cell density will converge towards zero if the oxygen concentration consistently falls below the survival threshold (42) Using the parameter values reported in Table 3 and the range of initial cell density reported in Table 6 (δ being an increasing function of n0) we find that this threshold ranges from cτ = 4.11%O2 for n0 = 1.00×107 cells/ml to cτ = 16.43%O2 for n0 = 4.00×108 cells/ml. The tissue being the only long-term source of oxygen with ctissue = 5%O2, the model will then predict no long-term cell survival for n0>3.53×107 cells/ml. This is in qualitative agreement with Fig 3A where only very low initial cell densities yield substantial cell survival. This behaviour is completely controlled by the parameter values, which were obtained by Coy et al. [12] using in vitro data corresponding to the first 24h. Hence, simulations beyond 24h should be considered as extrapolations and used as qualitative arguments.

With regards to the effect of cell seeding on VEGF secretion, Fig 3B displays the mean (black line), the range (light grey area) and the standard deviation (dark grey area) associated to VEGF concentration as a function of the initial cell density. The black dots represent the position of the initial cell density yielding the highest mean VEGF concentration and largest standard deviation at a given time point, i.e. and . In general, vSD represents an indirect measure of VEGF gradient in all directions. However, for the particular case of long and thin () cylinder surrounded by an impermeable sheath, radial gradients have only a small influence so that vSD mainly represents longitudinal gradients of VEGF. After 12h, the mean VEGF concentration globally exhibits variations similar to the mean cell density, with a maximum value for cells/ml, meaning that the highest initial cell density does not yield the largest VEGF secretion. High initial cell densities are associated with fast oxygen consumption, resulting in the cell secreting more VEGF at early time points. However, they are also associated with higher cell death (Fig 3A). In addition, effect of VEGF degradation, with characteristic time (43) cannot be neglected (Tv,K≈9.15h), implying that a VEGF molecule secreted early has likely been degraded after 12h, which is also detrimental to high initial cell density. For lower initial cell densities, on the other hand, the initial oxygen is consumed slower, so that the secretion of VEGF is delayed, hence the lower concentration after 12h. The balance between promoting VEGF secretion via quick oxygen consumption and maintaining a reasonable cell density to sustain meaningful VEGF secretion justifies the existence of a maximum in the mean VEGF concentration after 12h.

Fig 3B also shows that the lower bound of the VEGF concentration range is vanishingly small, which is a direct effect of the boundary condition (Eq 20). On the other hand, the upper bound globally exhibits the same variations as the mean VEGF concentration with a maximum for n0≈3.44×108cells/ml. This is expected since it represents the concentration of VEGF at the center of the construct, which behaves similarly to the mean VEGF concentration. The mean VEGF concentration is also generally closer to the upper bound of the range. The region perfused by oxygen from the tissue is relatively small (Eq 36) so that cell in the large central region are deprived of oxygen and secrete more VEGF. Given the slow axial diffusion, the concentration in VEGF in this region increases, and since it has a large volume, it contributes more to the mean VEGF concentration, which ends up being closer to the upper bound of the range.

Fig 3A finally shows that the standard deviation also globally follows the same variations as the mean and the upper bound of the VEGF concentration, with a maximum for cells/ml. We recall that the standard deviation serves as a global measure of the longitudinal VEGF concentration gradient, i.e., the difference between VEGF concentration at the centre (upper bound) and the extremities (lower bound). The latter being constant regardless of the initial cell density, the standard deviation is actually controlled by the upper bound, which explains its behaviour. The difference between and is due to the integration of axial non-uniformities in the VEGF concentration.

After 24h, Fig 3B shows that the mean VEGF concentration still exhibits the same variations with a maximum for cells/ml but has generally decreased except for n0≤2.5×107. Low initial cell densities are associated with slower oxygen consumption, sustaining longer cell proliferation (Fig 3A). This causes the VEGF secretion to peak later, explaining the increase in mean VEGF concentration. On the contrary, high initial cell densities are associated with an earlier peak in VEGF secretion. A large fraction of VEGF molecules secreted at that time have then likely been degraded after 24h, implying that the remaining concentration is the result of the much lower cell density (Fig 3A), explaining the decrease in the mean VEGF concentration. This balance between cell density, VEGF degradation and the peak in VEGF secretion also explains decrease.

Fig 3B also shows that the standard deviation exhibits the same behaviour as the mean VEGF concentration with cells/ml. We note that the difference between and has increased, which is the consequence of the continuing time integration of longitudinal non-uniformities in the VEGF concentration.

After 5 days, the VEGF mean, range and standard deviation have collapsed towards zero for most initial cell densities. At this time, only cells located close to the stumps remain alive. There, the oxygen level is above the hypoxic threshold so that the VEGF secretion is hindered. The combined effect of degradation and axial diffusion, which is no longer negligible (tv,L≈5.8 days), having eliminated the pre-existing VEGF concentration, explains the vanishingly small concentration.

Culture conditions.

In this section we quantify the impact of the initial oxygen (c0) and VEGF (v0) concentrations on cell survival and VEGF secretion. These conditions can be controlled experimentally to improve outcomes post-implantation. Fig 4A shows the initial cell density yielding the highest cell density (, dashed line), the highest VEGF concentration (, plain line) and the highest standard deviation in VEGF concentration (, dotted line) as functions of c0 after 24h. All three functions are decreasing, implying that lower initial cell densities yield better performances for high initial oxygen concentration. This is expected for since a higher initial oxygen concentration implies a longer consumption time (Eq 33) and therefore longer proliferation times. This situation favours lower initial cell densities after 24h (Fig 3A), explaining the behaviour of . These longer oxygen consumption and proliferation times also result in a delay in the peak of VEGF secretion, which also favours lower initial cell densities after 24h (Fig 3B), explaining the behaviour of and .

thumbnail
Fig 4. Increasing the initial oxygen concentration supports fewer seeded cells and results in an increased cell density and VEGF concentration and gradient after 24h.

A) initial cell density yielding the highest cell density (, dashed line), the highest VEGF concentration (, continuous line) and the largest standard deviation in the VEGF concentration (, dotted line) as functions of the initial concentration of oxygen i.e. c0. B) as a function of c0. C) (continuous line) and (dotted line) as functions of c0. D) (continuous line) and (dotted line) as functions of the initial concentration of VEGF i.e. v0. E) (continuous line) and (dotted line) as functions of v0.

https://doi.org/10.1371/journal.pcbi.1009142.g004

We note that so that there is no initial cell density that maximises both the mean VEGF concentration and the VEGF gradient. The differences between these two quantities are almost constant regardless of the initial oxygen concentration and is the result of the longitudinal non-uniformities in the VEGF concentration, which are not sensitive to the initial concentration of oxygen.

Fig 4B shows the maximum mean cell density, i.e., , as a function of the initial oxygen concentration in the construct. We can see that it increases with c0, which is expected since increasing the initial oxygen concentration is associated with an increase in the cell proliferation time resulting in an increase in the cell density after 24h.

Similarly, Fig 4C shows that the maximum mean VEGF concentration and maximum VEGF standard deviation, i.e., and , are also slightly increasing functions of c0. This is consistent with the increase in as a higher cell density, under low oxygen conditions, secretes more VEGF.

Finally, Fig 4D shows and and Fig 4E shows and , as functions of the initial VEGF concentration in the construct v0. We can see that v0 yields little to no effects on these quantities. This can be explained since most of the VEGF initially present has been degraded after 24h (Eq 43). In addition, ng/ml which is much higher than v0, explaining the constant behaviour of . The same reasoning applies for and recalling that vSD is mainly controlled by the concentration at the centre of the construct.

Surrounding tissue.

In this section we investigate the impact of the tissue oxygen concentration (ctissue) and tissue VEGF concentration (vtissue) upon cell survival and VEGF concentration. Both of these parameters will depend on the specific repair scenario and are challenging to define due to a paucity of direct measurements in the literature. Therefore, it is important to understand the impact of these parameters on model predictions. Fig 5A shows the initial cell density yielding the highest cell density (, dashed line), the highest VEGF concentration (, plain line) and the largest standard deviation in the VEGF concentration (, dotted line) as functions of ctissue after 24h. Similar to Fig 4A, is decreasing, so that lower initial cell densities perform best with a high tissue oxygen concentration. This can be explained since an increase in oxygen concentration, even when predominantly located close to the proximal and distal stumps, favours low initial cell densities (Fig 3A).

thumbnail
Fig 5. Increasing the surrounding oxygen tissue concentration supports fewer seeded cells and results in an increased cell density after 24h.

A) initial cell density yielding the highest cell density (, dashed line), the highest VEGF concentration (, continuous line) and the largest VEGF concentration standard deviation (, dotted line) as functions of the concentration of oxygen in the surrounding tissue, i.e., ctissue. B) as a function of ctissue. C) (continuous line) and (dotted line) as functions of ctissue. D) (continuous line) and (dotted line) as functions of the concentration of VEGF in the surrounding tissue i.e., vtissue on a logscale. E) (continuous line) and (dotted line) as functions of vtissue on a logscale.

https://doi.org/10.1371/journal.pcbi.1009142.g005

On the contrary, and are not sensitive to an increase in the tissue oxygen concentration. An increase in ctissue mainly reduces VEGF secretion by cells in the oxygen-perfused regions (Eq 36). In these regions, the VEGF concentration is already low, due to the combined effect of the tissue VEGF concentration (vtissue = 0) and longitudinal diffusion. Consequently, and vSD, which are mainly controlled by the VEGF concentration in the central region of the construct (Fig 3B), remain relatively unaffected by changes in ctissue.

Analogous to Fig 4B, Fig 5B shows that is an increasing function of ctissue, as cells close to proximal and distal stumps have access to a higher oxygen concentration and therefore proliferate more.

On the other hand, Fig 5C shows that and , are not sensitive to an increase in the tissue oxygen concentration. and vSD, and , are mainly controlled by the VEGF concentration in the central region of the construct, which is not significantly affected by an increase in ctissue.

Finally, Fig 5D and 5E show and respectively, as functions of the tissue concentration of VEGF. Similar to Fig 4D and 4E, none of them is sensitive to vtissue. This is because and are much higher than vtissue meaning that they are completely controlled by the VEGF secreted in the construct, also explaining the constant behaviour of and .

Porous sheath

In this section, we investigate the effect of sheath thickness, porosity and tortuosity upon cell survival, oxygen concentration and VEGF concentration. The goal is to use insights from the model to inform which sheath properties enable improved NRC performance after 1 day.

Fig 6 shows heatmaps of cell density, VEGF and oxygen concentrations over 24 hours. The cell density decreases over time despite cells having better access to oxygen, although this decrease is smaller than for the impermeable scenario. This decrease is due to the large cell death rate, which is not compensated by the improved access to oxygen. Contrarily to an impermeable sheath however, there are no discernible longitudinal heterogeneities in the cell density as ctissue is uniform so cells have similar access to oxygen regardless of their axial position. This holds only because the sheath has a large porosity (ϵ = 0.8) and a relatively small thickness (T = 0.25mm) so the radial diffusion time of oxygen in the construct

thumbnail
Fig 6. Cell-solute interactions in a permeable nerve repair construct results in a quick decrease in oxygen concentration and a low uniform VEGF concentration.

Heatmaps of cell density, oxygen and VEGF concentration over 24 hours for an initial cell density n0 = 1.78×106 cells/ml and a porosity ϵ = 0.8. The porous sheath structure is approximated by an array of overlapping spheres (see Table 1) with a thickness T = 0.25mm and a porosity ϵ = 0.8.

https://doi.org/10.1371/journal.pcbi.1009142.g006

(44) is much shorter (tc,R≈0.05h) than its axial counterpart (tc,L≈34.5h≈1.44 days).

The VEGF concentration evolves similarly, although much faster, to the case of an impermeable sheath (Fig 2). The concentration is much lower with no discernible axial gradients. This is due to the combined effect of zero tissue VEGF concentration and the radial diffusion time of VEGF in the construct (45) being much shorter (tv,R≈0.31h) than its longitudinal counterpart (tv,L≈139h≈5.8 days), and the secretion time of VEGF (tv,P≈3.63h). This equilibrium between radial diffusion and VEGF secretion explains why VEGF concentration rises to such low values at early times, and the difference between radial and longitudinal diffusion times explains why there are no longitudinal gradients. The long-term decrease of VEGF is due to the cell density decrease combined with the oxygen concentration in the tissue being higher than the hypoxic threshold (ctissue>ch).

Finally, Fig 6 shows that oxygen concentration relaxes from c0 to ctissue in under 0.5h. Such a quick relaxation is due to the combination of the cell density being the highest at early times, which implies stronger oxygen consumption, with the radial diffusion time of oxygen in the construct being small (tv,R≈0.05h).

Fig 6 focuses on a single initial cell density and a single porosity and does not quantitatively inform the best cell seeding strategy. By comparison, Fig 7 explores the impact of seeded cell density for different sheath porosities. Consistent with Fig 3A, the oxygen concentration decreases faster for higher initial cell densities when the sheath is impermeable. We also note that when the sheath is porous the concentration in oxygen remains globally constant regardless of the initial cell density and is close to ctissue in under 0.5h, even for ϵ = 0.1, which is associated with a smaller radial diffusive flux (Eq 15). This is a direct consequence of the radial diffusion time of oxygen (tc,R<0.1h) being much smaller than the consumption time of oxygen (0.96h<tc,M<38h), so that the radial diffusive flux is not limiting.

thumbnail
Fig 7. Increasing the sheath porosity of a permeable construct results in a higher cell density and oxygen concentration and a much lower VEGF concentration and gradient after 24h.

Mean cell density , mean oxygen concentration , mean VEGF concentration and standard deviation vSD, as functions of the initial cell density n0, for different sheath porosity ϵ, at different time points. The porous sheath structure is approximated by an array of overlapping sphere (see Table 1) with a thickness T = 0.25mm. We note that ϵ = 0 corresponds to an impermeable sheath (red lines).

https://doi.org/10.1371/journal.pcbi.1009142.g007

With regard to cell survival, we see that for t<1h the cell density associated with the impermeable sheath is slightly higher than the cell densities associated with a porous sheath, regardless of n0. This is due to the quick loss of oxygen by radial diffusion when the sheath is porous. Even after 12h, for n0≤4×107cells/ml, the impermeable sheath still performs slightly better than the porous sheath despite exhibiting lower oxygen concentration. That is because the oxygen concentration only affects the cell proliferation, i.e. the derivative of the cell density, and not the cell density directly. This introduces a delay between the moment when the oxygen concentration in the impermeable case becomes smaller than in the porous case, and the moment when the cell density effectively becomes lower than in the porous case. After 24h, the oxygen concentration in the impermeable case has been lower long enough so that the associated mean cell density is also lower than the porous case, regardless of n0.

After 12h and 24h both and are higher when the sheath is porous ( cells/ml and cells/ml after 24h for ϵ≥0.1) as the cells have access to oxygen supplied by the tissue. Finally, we note that the cell density has the same behaviour regardless of the porosity, which aligns with the mean oxygen concentration also being insensitive to porosity.

With regards to VEGF secretion, a porous sheath dramatically decreases VEGF concentration after 24h compared to an impermeable sheath. However, after 1h both the mean and standard deviation of the VEGF concentration are higher for a porous sheath with ϵ = 0.1 than for an impermeable sheath, except for cell densities above n0>3.30×108cells/ml. This is because 1) the oxygen concentration is lower when the sheath is porous except for n0>3.30×108cells/ml, so that cells will secrete more VEGF, 2) ϵ = 0.1 implies a small radial diffusive flux between the collagen cylinder and the annular sheath (Eq 16) and 3) the radial diffusion time is long enough (tv,R = 0.63h) so that VEGF concentration can build inside the construct. In line with this observation, increasing the porosity leads to a significant decrease in the mean VEGF concentration and standard deviation.

Sheath material.

In this section, we investigate the effect of sheath tortuosity, porosity and thickness on cell survival and VEGF secretion. Fig 8A shows, for each porous structure introduced in Table 1, and the associated maximum cell density as functions of the sheath porosity and thickness after 24h. We see that for porosity in the range 0.3≤ϵ≤1 all three porous structures yield similar performances. For smaller porosity (ϵ<0.3), however, we see that the cylindrical structure differs significantly from the two others. That is because the associated tortuosity model has been asymptotically derived for large porosity and exhibits a singularity at ϵ = 0.25. We also note that increasing the porosity results in a sharp increase of and for the spherical and tubular structures that plateaus for ϵ≥0.3 because cells have better access to oxygen via radial diffusion. Small porosities are associated with a small radial flux of oxygen ( in Eq 15) so that increasing ϵ relaxes such constraint.

thumbnail
Fig 8. Increasing sheath porosity increases the maximum cell density and decreases the VEGF concentration and gradient field after 24h, regardless of the underlying archetypal porous structure.

Increasing sheath thickness has opposite effects. A) and the associated initial cell density and B) and the associated initial cell density as functions of sheath porosity ϵ and sheath thickness T for the different tortuosity expressions listed in Table 1.

https://doi.org/10.1371/journal.pcbi.1009142.g008

Finally, we see that increasing the sheath thickness decreases and noticeably for ϵ≤0.3. An increase in sheath thickness implies an increase in radial diffusion time of oxygen (e.g. tc,R≈0.05h for T = 0.25mm and tc,R≈0.5h for T = 1.5mm) so that cells have an impaired access to oxygen, especially for small porosities, which are already associated with a limited flux of oxygen ( in Eq 15). This explains the decrease in and .

Fig 8B shows and as functions of the sheath porosity and thickness after 24h. Results with regard to the mean VEGF concentration were not reported here as they were very similar. Both and are globally not sensitive to the underlying porous structure. Contrarily to cell density, however, increasing sheath porosity decreases and . This is due to the combined effect of the relatively short radial diffusion time and the increase in the flux of VEGF at the interface between collagen and sheath.

Finally, increasing sheath thickness also increases and , especially for small sheath porosity. This is expected as a thicker sheath implies a longer radial diffusion time for VEGF (e.g., tv,R≈0.31h for T = 0.25mm and tv,R≈5.22h for T = 1.5mm, assuming ϵ = 0.8 and a spherical porous structure) so that molecules can accumulate in the construct. Increasing sheath thickness also impairs access to oxygen in the main body of the NRC, so that cells rely on the initial oxygen concentration to survive. Higher initial cell densities tend to produce more VEGF due to quick oxygen consumption, hence the increase in .

Partially porous sheath.

In this section, we investigate the impact of a sheath that is porous at the extremities and impermeable at the centre. This is motivated by the idea that a partially porous sheath could serve as a compromise between the impermeable case that favours VEGF secretion, and the porous case that favours cell survival.

Fig 9A shows heatmaps of cell density and VEGF concentration after 24h, for two partially porous sheaths (ϵ = 0.8 and T = 0.25mm) with porous length Lp = 1mm and Lp = 3mm respectively. We can see that the porous regions exhibit higher cell density and lower VEGF concentration compared to the impermeable region, which is consistent with results reported in Fig 7, and is due to short radial diffusion times and the cells having improved access to oxygen. By comparison, the central region has only access to oxygen through longitudinal diffusion, which is still significantly slower (tc,L≈1.1 days for Lp = 1mm and tc,L≈12.5h for Lp = 3mm substituting L by L−2Lp in Eq 32) than the consumption of oxygen (e.g. for n0 = 1.78×108 cells/ml, tc,M≈1.95h), so that cells in this region mostly rely on the initial oxygen supply. This initial oxygen reserve is consumed within a few hours, leading to lower oxygen concentration and lower cell density after 24h. By comparison, the slow longitudinal diffusion timescale for VEGF in the central region (tv,L≈4.3 days for Lp = 1mm and tv,L≈2 days for Lp = 3mm substituting L by L−2Lp in Eq 34) together with the low oxygen concentration there gives rise to a higher VEGF concentration.

We still note that increasing Lp results in an increase in cell density and a decrease in VEGF concentration in the central region, mainly due to the reduction in oxygen and VEGF longitudinal diffusion times.

thumbnail
Fig 9. A partially porous sheath with a moderate length Lp limits the decrease of VEGF concentration while still increasing cell density and oxygen concentration after 24h.

A) Heatmaps of cell density and VEGF concentration after 1 day, for two different porous sheath zone lengths Lp, for an initial cell density n0 = 1.78×106 cells/ml. B) Mean cell density , mean oxygen concentration and mean VEGF concentration and standard deviation vSD, as functions of the initial cell density n0, for different Lp values over 1 day. In A) and B) the porous sheath structure is approximated by an array of overlapping spheres (see Table 1) with a thickness T = 0.25mm and a porosity ϵ = 0.8. We note that Lp = 0 (red lines) corresponds to an impermeable sheath.

https://doi.org/10.1371/journal.pcbi.1009142.g009

Fig 9B extends these results and shows the mean cell density, VEGF and oxygen concentration, as well as the standard deviation in the VEGF concentration, as functions of the initial cell density for different length of the porous sheath zones (Lp). After 0.5h, the mean oxygen concentration decreases with Lp regardless the initial cell density. This is because 1) oxygen concentration in the porous regions quickly decreases to reach ctissue and 2) increasing Lp decreases the longitudinal diffusion time in the impermeable region (substituting L by L−2Lp in Eq 32), resulting in the observed decrease in mean oxygen concentration. After 1h, however, we see that the mean oxygen concentration starts to increase with the porous length for high initial cell densities (n0>3.30×108cells/ml). This is a consequence of the short oxygen consumption time, which causes the oxygen concentration in the impermeable region to decrease below the concentration in the porous region. Increasing Lp hence increases the mean oxygen concentration by reducing both the size of the impermeable region and the longitudinal diffusion time. This also explains why after 24h increasing Lp increases the mean oxygen concentration regardless of the initial cell density.

Fig 9B also shows that after 24h the mean cell density increases with Lp, independent of the initial cell density. This is consistent with the mean oxygen concentration increasing with Lp, and consistent with Fig 7 for a uniformly porous sheath. Still, we note that for the first few hours the variations are slightly more complex due to the high initial oxygen concentration, which tends to favour low initial cell densities when the sheath is impermeable.

Fig 9B also shows that increasing Lp decreases the mean VEGF concentration regardless of time or initial cell density. This is the result of 1) the low VEGF concentration in the porous region, 2) the increase in oxygen concentration in the impermeable region, which hinders VEGF secretion and 3) the decrease of the VEGF longitudinal time of diffusion in the impermeable region.

For vSD, on the other hand, Fig 9B shows that variations are more complex and non-monotonic. Generally, we see that increasing Lp first increases the standard deviation to a maximum value and then decreases it. When the sheath is completely impermeable (Lp = 0), the mean VEGF concentration is much closer to the maximum of concentration (Fig 3B), which means that the VEGF concentration it is generally distributed close to its mean. Increasing Lp then decreases the mean VEGF concentration, but not so much the maximum VEGF concentration at the centre of the impermeable region, mainly due to slow axial diffusion and limited access to oxygen. This causes the VEGF concentration distribution to spread more, increasing vSD. Increasing Lp further results in the mean VEGF concentration being closer to the concentration in porous regions and a decrease in the maximum VEGF concentration due to better access to oxygen and shorter longitudinal diffusion in the impermeable region, effectively decreasing vSD.

Non-uniform cell seeding strategies

In this section, we investigate the effect of different non-uniform initial cell distributions upon cell survival and VEGF secretion. The premise is that supporting denser populations of cells in specific regions could induce a differential in VEGF secretion across the construct, thereby enhancing VEGF gradients which act as a stimulus for vascular regrowth.

To create non-uniformities in the cell distribution we consider the following expression for the initial cell density (46) where n0,tot denotes the initial total number of cells seeded and ζ the repartition factor controlling the magnitude of the non-uniformity across 3 equally spaced regions. A small repartition factor (ζ<1) concentrates the cells in the regions adjacent to the stumps while a larger one (ζ>1) concentrates the cells in the central region of the construct (Fig 10A). Although we focus on dividing the NRC into 3 zones, the approach can be extended. We note that each value of n0,tot is related to a uniform cell seeding density n0 by assuming ζ = 1. For instance, n0≈170, 102 or 34×106cells/ml yield n0,tot = 5.00, 3.00 or 1.00×105 cells respectively. Considering the above square function for the initial cell density has two advantages. First, its simplicity will facilitate ease of manufacture. Second, integrating it over the whole construct consistently yields the initial total number of cells, regardless of the repartition factor.

thumbnail
Fig 10. Overconcentrating cells in the centre of the construct leads to lower cell density and VEGF concentration after 24h.

A) Heatmaps of initial cell density for different repartition factor ζ for a constant initial total number of cells n0,total = 2×105 cells. The case ζ = 1 corresponds to a uniform distribution. B) Heatmaps of cell density (left) and VEGF concentration (right) for impermeable (top) and porous (bottom) sheaths after 1 day. In each panel, simulations are run for two different values of ζ and n0,total. When porous, the sheath structure is approximated by an array of overlapping spheres (see Table 1) with a thickness T = 0.25mm and a porosity ϵ = 0.8.

https://doi.org/10.1371/journal.pcbi.1009142.g010

Fig 10B shows heatmaps of cell densities and VEGF concentrations after 24h, for different initial total number of cells, for a uniformly porous and an impermeable sheath, and for uniform (ζ = 1) and centre-populated (ζ = 3) cases. On the top left panel, the sheath is impermeable, and we see that when ζ = 3 and ntot,0 = 1.00×105cells, the cell density in the central region is higher than in the stump regions. Comparing it with the uniform case, we also see that the cell density in the central region is higher and that the cell densities in the stump regions are lower.

To help interpret such results we compute, for the central and stump regions, the associated initial cell densities. We then compare them to , obtained in Fig 3A for a uniform seeding, with the idea that each region can approximately be treated independently so that can serve as a proxy to represent the balance between effects of initial cell density, cell proliferation and cell death. The cell density in the central region for the case ζ = 3 and ntot,0 = 1.00×105cells is n0≈6.10×107cells/ml, which is closer to cells/ml than the initial cell density in the regions adjacent to the stumps (n0≈2.00×107cells/ml) or in the uniform case (n0≈3.40×107cells/ml), explaining the higher cell density in the central region.

By comparison, we see that when ζ = 3 and ntot,0 = 5.00×105cells, the cell density in the central region is much lower. This time n0≈3.05×108cells/ml in this region, which is relatively far from , and hence is associated with low cell density. On the contrary n0≈1.01×108cells/ml in the stump regions, which is significantly closer to explaining the higher cell density, even compared to ζ = 1 (n0≈1.70×108cells/ml).

The top right panel shows that the corresponding VEGF concentration distributions remain smooth despite discontinuity in the cell density thanks to longitudinal diffusion. Still, we see that when ntot,0 = 1.00×105cells the VEGF concentration is higher when the cells are concentrated in the central region of the construct. On the contrary, when ntot,0 = 5.00×105cells the VEGF concentration is higher when the cells are seeded uniformly, although we note that the VEGF concentration still increases between ntot,0 = 1.00×105cells and ntot,0 = 5.00×105cells for both cases. Similar to the cell density, such results can be interpreted by comparing to the value of determined in Fig 3B. For instance, ntot,0 = 5.00×105cells with ζ = 1 yields n0≈1.70×108cells/ml which is closer to cells/ml than any other configuration, explaining the higher VEGF concentration.

The bottom left panel shows that the cell density is globally higher and follows the same pattern when the sheath is porous, regardless of the initial cell density or seeding strategy. This is mainly due to the cells having improved access to oxygen via radial diffusion.

Finally, the bottom right panel shows that when the sheath is porous, VEGF concentrations are very low, regardless of the cell-seeding strategy. This is primarily due to the VEGF radial diffusion time being smaller than the VEGF secretion time (tv,R≈0.31h versus tv,P≥1.65h for the 4 cases considered here) and the cells having improved access to oxygen from the tissue.

Fig 11A extends these results and shows ntot after 12h and 24h, as a function of ζ, for different and for a porous and an impermeable sheath.

thumbnail
Fig 11. Seeding cells in the stump regions increases cell density and seeding cells in the central region increases VEGF concentration after 24h for an impermeable sheath.

For a porous sheath, uniform seeding leads to maximal cell density. A) Total cell number and B) mean, standard deviation (dark grey areas) and ranges (light grey areas) of VEGF concentration, as functions of the repartition factor ζ and the initial total number of cells n0,total, for an impermeable (continuous line) and a porous (dashed line) sheath, after 12 and 24 hours. Black dots indicate the repartition factor ζ yielding the highest total number of cells, i.e. the highest VEGF mean concentration, i.e. and the highest standard deviation of VEGF concentration, i.e. . When porous, the sheath structure is approximated by an array of overlapping spheres (see Table 1) with a thickness T = 0.25mm and a porosity ϵ = 0.8.

https://doi.org/10.1371/journal.pcbi.1009142.g011

When the sheath is porous, we see that the total number of cells is non monotonic and maximal for repartition factors close to ζ = 1, so that concentrating cells either in the proximal and distal stumps regions or in the central region of the NRC is not beneficial. In this situation the oxygen concentration in the construct is approximately uniform due to quick radial diffusion so that concentrating cells in a specific region only increases the associated initial cell density, which in turn increases the cell death rate. Therefore, uniform seeding remains the best strategy as it allows maximum spreading.

When the sheath is impermeable, on the other hand, we see that concentrating cells in the stump regions is beneficial for small . As the initial cell number increases, however, we see that increasing ζ, i.e., the number of cells in the central region, increases ntot. This is mainly the result of the high initial oxygen concentration that can sustain cell proliferation even when cells are concentrated in a given region. However, seeding too many cells in a specific region, again results in a smaller total number of cells after 24h.

Fig 11A also shows that after 12h, the total number of viable cells is larger when the sheath is impermeable, i.e. for cells and for cells, which are associated with relatively low initial cell densities (n0≤3.10×107cells/ml for cells and n0≤6.20×107cells/ml for cells when ζ = 3). This is again an effect of the high initial oxygen concentration which enhances cell proliferation when the sheath is impermeable. In all other instances, in particular after 24h, the porous sheath exhibits larger total cell numbers because the cells have an improved access to oxygen.

Fig 11A shows the existence of an initial total cell number beyond which the total number of cells decrease after 24h. The precise position of such a maximum slightly depends on the repartition factor and whether the sheath is porous or impermeable, but we see that it is located around cells (equivalent to n0 = 1.02×108cells/ml). This aligns with values of reported in Fig 3A (8.80×107cells/ml) for an impermeable sheath and in Fig 7 (1.10×108cells/ml) for a porous sheath.

With regard to VEGF, Fig 11B shows the mean, the standard deviation, and the range of VEGF concentration, as functions of the repartition factor ζ, for different initial total number of cells and for an impermeable and a porous sheath. Contrary to ntot, we see that seeding the cells in the central region increases both the mean VEGF concentration and the standard deviation when the sheath is impermeable. The central region is associated with the longest longitudinal diffusion time so that seeding the cells preferentially in this region allows the VEGF concentration to rise to higher values. Still, we note that when the initial total cell number increases it becomes preferable to seed the cells more uniformly. This is the result of the balance between oxygen consumption, cell death and VEGF secretion. When the initial total number of cells is relatively small, the initial concentration of oxygen sustains the cell population longer even if the cells are concentrated in the central region. As the initial total number of cells increases, the oxygen is consumed faster and the cell death rate increases, especially in the central region, which leads to lower cell density, ultimately reducing VEGF secretion. Spreading the cells then increases the consumption time of oxygen while decreasing the cell death rate, sustaining a higher VEGF production.

By comparison, increasing generally results in an increase of the VEGF production after 24h. This is in agreement with Fig 3B which shows that and are significantly larger than and correspond to greater than the values explored in Fig 11B.

After 12h, we see that the maximum of the VEGF concentration range exhibits two maxima (ζ = 0 and ζ = 3) and one minimum (ζ≈1). These two peaks are the results of the cells being concentrated in specific regions. The minimum therefore corresponds to the transition between these two modes of VEGF secretion. We note that after 24h, for cells, this behaviour is reverted, with two minima around ζ = 0 and ζ = 3 and a maximum value around ζ = 1. In this scenario, the oxygen initially present in the seeded region is quickly consumed which results in a quick decay of the associated cell population. Finally, we note that the VEGF standard deviation also follows such behaviour. This is expected as we have shown that vSD is controlled by the maximum VEGF concentration (Fig 3B).

Consistent with previous analysis, we see that when the sheath is porous, the mean, standard deviation and range of VEGF concentration are all much lower regardless of the repartition coefficient. This is due to the quick radial diffusion of VEGF and the cells having better access to oxygen, which prevents the VEGF concentration from increasing even when the cells are concentrated in the central region.

Altogether Fig 11 shows that non-uniform seeding does not change dramatically the behaviours described in Figs 3 and 7 for uniform cell seeding.

Discussion

We devised and simulated a cell-solute model with the goal to comparatively evaluate the short-term performance of a range of NRC designs and to link these outcomes to the underlying transport mechanisms. This approach, using mathematical models informed by experimental data, may identify design strategies for NRCs as well as improve our understanding of the cellular and chemical factors which determine outcomes. Given the challenges and cost of evaluating the full design parameter space using a purely experimental approach, computational modelling may be valuable in refining design choices and informing future experimental work.

Our primary focus was to determine scenarios which promote high cell viability, as well as VEGF concentration and gradients, given that these are key determinants of vascular growth (required to sustain an implanted construct in the longer term). We focused on critical design features which may be controlled when fabricating engineered NRCs: seeded cell density and spatial distribution, and the permeability of the enclosing NRC sheath to molecular factors. The cell-solute model was parameterised against available literature data, including dedicated in vitro experiments reported in [12] on a therapeutically-relevant cell type, dADSCs.

Model predictions showed that seeding more cells in a NRC does not necessarily lead to higher cell densities after 24h due to the balance between oxygen transport, cell proliferation and cell death. This result runs contrary to the reasoning generally used in experimental studies, which either neglects to consider the effect of seeding cell densities and distributions entirely or assumes that seeding more cells will necessarily result in a greater number of viable cells over time. In the case of an impermeable sheath, we found that the initial cell density yielding the highest cell density after 24h is cells/ml, which corresponded to less than half the maximal initial cell density tested (Fig 3A). Interestingly, this fell close to the value of 8.00×107cells/ml found by Mosahebi et al. [35] to produce on average the greatest length of axonal regeneration after 3 weeks in vivo in a rat sciatic nerve injury model. Increasing the cell density beyond this showed a decrease in regeneration in their experimental study. Guénard et al. [36] produced similar experimental results, also in rat sciatic nerve models, concluding that 1.20×108cells/ml produced the greatest number of myelinated axons after 3 weeks, although in this case the luminal diameter was altered along with the cell seeding density.

We also showed that when the sheath was impermeable was noticeably lower than , so that it was not possible to maximise both the VEGF production and the cell density at the same time (Fig 3). We note however that and were still relatively close to each other so that a seeding strategy maximising the mean VEGF concentration also enabled steep VEGF gradients after 24h (Fig 3B).

We further pointed out that for cells/ml the associated cell density was cells/ml, which yielded a relatively small survival rate after 24h (Fig 3A). Such a cell density is relatively low and may not provide the required support for nerve regeneration over the longer term (e.g. angiogenesis, neurite growth), or sustain the long-term secretion of VEGF.

We showed that seeding the cell preferentially in the stump regions could help increase the viable cell number after 24h (Fig 11A) and that seeding cells in the centre of the construct on the contrary increased the VEGF concentration and gradients (Fig 11B). We note, however, that in each scenario the optimal repartition factor ζ remained close to unity, so that non-uniform seeding strategies did not considerably improve VEGF production or cell density after 24h.

By comparison, increasing the ambient construct oxygen concentration at the outset (via the culture conditions) appeared to be more beneficial, as it resulted in an increase of and , and the decrease of (Fig 4A-C). Essentially a high culture oxygen concentration allowed fewer cells to be seeded, which then proliferated well and the viable cell density after 24h was systematically improved along with the secretion of VEGF.

Increasing the initial concentration of VEGF yielded no significative changes to the VEGF field after 24h, mainly due to the range explored being small compared to VEGF concentration originating from cell production. This suggests alternative approaches would be needed to boost VEGF levels in NRCs, for example controlled-release VEGF via biomaterials or nanoparticles. In these cases, the concentration of delivered VEGF should be higher than the range simulated here and could be spatially distributed to augment both VEGF concentration and gradients.

When the sheath was porous, however, the initial concentration of oxygen had little influence on the results due to quick radial diffusion, so that the problem was controlled by the ambient oxygen levels in the surrounding tissue (Fig 7). On one hand, this increased viable cell density after 24h yielding better cell survival (), because cells had improved access to oxygen. On the other hand, this dramatically decreased VEGF concentration and gradients, even when the sheath porosity was relatively small mainly due to quick radial diffusion and the zero VEGF concentration in the tissue, but also because the surrounding tissue concentration of oxygen was above the threshold for VEGF release. Although an option may be to increase the sheath thickness while keeping the sheath porosity relatively low (Fig 8B), this also impaired oxygen availability and resulted in a decrease of viable cell density after 24h (Fig 8A).

A partially porous sheath with an intermediate value for the length of the porous zones (e.g. 2mm<Lp<4mm) appears to address the compromise between sustaining cell viability and stimulating VEGF concentration and gradients (Fig 9B). However, this raises another consideration. Endothelial cells forming the tip of vessel sprouts invading the NRC during angiogenesis will sense VEGF thanks to filopodia extensions that can extend up to 100 microns [74]. A partially porous sheath with a porous zone length of 2-4mm would create a region depleted of VEGF much larger than the range of filopodia, which might hinder endothelial cell invasion and angiogenesis. The balance of these factors should be explored in future experimental work.

Overall, predictions of the computational model indicate that a porous sheath has both advantages and drawbacks. This ambiguity is also present in experimental studies; some report increased neuronal regeneration [7578] and offer the hypothesis that this could be a result of increased oxygen diffusion and infiltration of supportive cells, whereas others report poorer results when using porous materials [62,79] and suggest this is caused by loss of growth factors and infiltration of host immune cells. We point out however that the model developed in this work is limited to archetypal porous structures and offers only a simplified description of the initial steps of the nerve repair process. Further analysis, combining model predictions with experimental investigation, is required to explore this further.

A way to further assess the benefits of a porous sheath could be to extend simulations beyond the first 24h. The model should then include a description of angiogenic and neurogenic processes taking place on the longer timescale. Similarly, changes in tissue concentration in oxygen and growth factors are likely to happen beyond 24h and should also be included in the model given the importance of the tissue oxygen concentration on determining cell outcomes in the NRC (Fig 5A and 5B). Additionally, an empirical characterisation of the porous structure of the NRC sheath used for EngNT NRC would improve the model physical accuracy. Such an extended model will comprise an array of new parameters, some of which being very hard to estimate as they are not directly measured in vivo (e.g. evolving oxygen concentration in the tissue, endothelial sprouting growth rate) and will require dedicated experiments to narrow down their range.

Finally, an interesting compromise may be a porous sheath with selective properties that take advantage of the difference in molecular weight between VEGF (46kDa for a VEGF dimer [80]) and oxygen (3.2×10−2kDa), with pores smaller than VEGF molecules, but larger than oxygen. Such a porous sheath could then prevent loss of VEGF by radial diffusion whilst enabling diffusive transport of oxygen.

In conclusion, the simulation results presented here can be used to discard NRC designs that lead to sub-optimal cell survival and VEGF secretion within the first 24 hours post-implantation. In particular, results indicate that seeding cells beyond a given threshold is detrimental regardless of the situation, that non-uniform seeding is generally not significantly beneficial and that a porous sheath prevents the build-up of a significant VEGF concentration in the NRC. These predictions can now be used to propose new candidate NRC designs for in vivo testing. Finally, this study constitutes a first step in laying down the theoretical framework to accelerate the development of treatment strategies for peripheral nerve injuries.

References

  1. 1. Chen S-l, Chen Z-G, Dai H-l, Ding J-x, Guo J-s, Han N, et al. Repair, protection and regeneration of peripheral nerve injury. Neural regeneration research. 2015;10(11):1777. pmid:26807113
  2. 2. Grinsell D, Keating CP. Peripheral Nerve Reconstruction after Injury: A Review of Clinical and Experimental Therapies. BioMed Research International. 2014;2014:698256. pmid:25276813
  3. 3. Kallio PK, VastamÄKi M. An Analysis of the Results of Late Reconstruction of 132 Median Nerves. Journal of Hand Surgery. 1993;18(1):97–105. pmid:8436871
  4. 4. Lee SK, Wolfe SW. Peripheral Nerve Injury and Repair. JAAOS—Journal of the American Academy of Orthopaedic Surgeons. 2000;8(4). pmid:10951113
  5. 5. Vastamäki M, Kallio PK, Solonen KA. The results of secondary microsurgical repair of ulnar nerve injury. The Journal of Hand Surgery: British & European Volume. 1993;18(3):323–6. pmid:8393906
  6. 6. Wood MB. Peroneal Nerve Repair: Surgical Results. Clinical Orthopaedics and Related Research®. 1991;267. pmid:2044280
  7. 7. Carvalho CR, Oliveira JM, Reis RL. Modern Trends for Peripheral Nerve Repair and Regeneration: Beyond the Hollow Nerve Guidance Conduit. Frontiers in Bioengineering and Biotechnology. 2019;7(337). pmid:31824934
  8. 8. Wilcox M, Gregory H, Powell R, Quick TJ, Phillips JB. Strategies for Peripheral Nerve Repair. Current Tissue Microenvironment Reports. 2020;1(2):49–59. pmid:33381765
  9. 9. Melanie Georgiou SCJB Heather A. Davies, Loughlin Alison J., Golding Jonathan P., Phillips James B. Engineered neural tissue for peripheral nerve repair. Biomaterials. 2013;34(30):7335–43. pmid:23834895
  10. 10. Lackington WA, Ryan AJ, O’Brien FJ. Advances in Nerve Guidance Conduit-Based Therapeutics for Peripheral Nerve Repair. ACS Biomaterials Science & Engineering. 2017;3(7):1221–35. pmid:33440511
  11. 11. Georgiou M, Golding JP, Loughlin AJ, Kingham PJ, Phillips JB. Engineered neural tissue with aligned, differentiated adipose-derived stem cells promotes peripheral nerve regeneration across a critical sized defect in rat sciatic nerve. Biomaterials. 2015;37:242–51. pmid:25453954
  12. 12. Coy R, Al Badri G, Kayal C, Kingham P, Phillips J, Shipley R. Combining in silico and in vitro models to inform cell seeding strategies in tissue engineering. Journal of The Royal Society Interface. 2020. pmid:32208821
  13. 13. O’Rourke C, Day AGE, Murray-Dunning C, Thanabalasundaram L, Cowan J, Stevanato L, et al. An allogeneic ‘off the shelf’ therapeutic strategy for peripheral nerve tissue engineering using clinical grade human neural stem cells. Scientific Reports. 2018;8(1):2951. pmid:29440680
  14. 14. Carvalho CR, Reis RL, Oliveira JM. Fundamentals and Current Strategies for Peripheral Nerve Repair and Regeneration. In: Chun HJ, Reis RL, Motta A, Khang G, editors. Bioinspired Biomaterials: Advances in Tissue Engineering and Regenerative Medicine. Singapore: Springer Singapore; 2020. p. 173–201. pmid:32602098
  15. 15. Waters SL, Schumacher LJ, El Haj AJ. Regenerative medicine meets mathematical modelling: developing symbiotic relationships. npj Regenerative Medicine. 2021;6(1):24. pmid:33846347
  16. 16. Coy RH, Evans OR, Phillips JB, Shipley RJ. An integrated theoretical-experimental approach to accelerate translational tissue engineering. Journal of Tissue Engineering and Regenerative Medicine. 2018;12(1):e53–e9. pmid:27792286
  17. 17. An G, Mi Q, Dutta-Moscato J, Vodovotz Y. Agent-based models in translational systems biology. WIREs Systems Biology and Medicine. 2009;1(2):159–71. pmid:20835989
  18. 18. Metzcar J, Wang Y, Heiland R, Macklin P. A Review of Cell-Based Computational Modeling in Cancer Biology. JCO Clinical Cancer Informatics. 2019;(3):1–13. pmid:30715927
  19. 19. Van Liedekerke P, Palm MM, Jagiella N, Drasdo D. Simulating tissue mechanics with agent-based models: concepts, perspectives and some novel results. Computational Particle Mechanics. 2015;2(4):401–44.
  20. 20. Truesdell C, Noll W. The Non-Linear Field Theories of Mechanics. S.S. A, editor. Berlin, Heidelberg: Springer; 2004.
  21. 21. Whitaker S. The method of volume averaging: Springer Science & Business Media; 1981.
  22. 22. O’Dea R, Byrne H, Waters S. Continuum Modelling of In Vitro Tissue Engineering: A Review. In: Geris L, editor. Computational Modeling in Tissue Engineering. Berlin, Heidelberg: Springer Berlin Heidelberg; 2013. p. 229–66.
  23. 23. Lewis MC, MacArthur BD, Malda J, Pettet G, Please CP. Heterogeneous proliferation within engineered cartilaginous tissue: the role of oxygen tension. Biotechnology and Bioengineering. 2005;91(5):607–15. pmid:16025534
  24. 24. Malda J, Rouwkema J, Martens DE, le Comte EP, Kooy FK, Tramper J, et al. Oxygen gradients in tissue-engineered Pegt/Pbt cartilaginous constructs: Measurement and modeling. Biotechnology and Bioengineering. 2004;86(1):9–18. pmid:15007836
  25. 25. Obradovic B, Meldon JH, Freed LE, Vunjak-Novakovic G. Glycosaminoglycan deposition in engineered cartilage: Experiments and mathematical model. AIChE Journal. 2000;46(9):1860–71.
  26. 26. Shipley RJ, Waters SL. Fluid and mass transport modelling to drive the design of cell-packed hollow fibre bioreactors for tissue engineering applications. Mathematical Medicine and Biology: A Journal of the IMA. 2011;29(4):329–59. pmid:22076984
  27. 27. Pohlmeyer JV, Waters SL, Cummings LJ. Mathematical Model of Growth Factor Driven Haptotaxis and Proliferation in a Tissue Engineering Scaffold. Bulletin of Mathematical Biology. 2013;75(3):393–427. pmid:23358798
  28. 28. Lemon G, King JR, Byrne HM, Jensen OE, Shakesheff KM. Mathematical modelling of engineered tissue growth using a multiphase porous flow mixture theory. Journal of Mathematical Biology. 2006;52(5):571–94. pmid:16463188
  29. 29. O’Dea RD, Waters SL, Byrne HM. A multiphase model for tissue construct growth in a perfusion bioreactor. Mathematical Medicine and Biology: A Journal of the IMA. 2009;27(2):95–127. pmid:19805485
  30. 30. Lemon G, King JR. Multiphase modelling of cell behaviour on artificial scaffolds: effects of nutrient depletion and spatially nonuniform porosity. Mathematical Medicine and Biology: A Journal of the IMA. 2007;24(1):57–83. pmid:17018570
  31. 31. Ochoa JA, Stroeve P, Whitaker S. Diffusion and reaction in cellular media. Chemical Engineering Science. 1986;41(12):2999–3013. https://doi.org/10.1016/0009-2509(86)85036-9.
  32. 32. Wood BD, Whitaker S. Diffusion and reaction in biofilms. Chemical Engineering Science. 1998;53(3):397–425. https://doi.org/10.1016/S0009-2509(97)00319-9.
  33. 33. Chung CA, Yang CW, Chen CW. Analysis of cell growth and diffusion in a scaffold for cartilage tissue engineering. Biotechnology and Bioengineering. 2006;94(6):1138–46. pmid:16586509
  34. 34. Chung CA, Chen CW, Chen CP, Tseng CS. Enhancement of cell growth in tissue-engineering constructs under direct perfusion: Modeling and simulation. Biotechnology and Bioengineering. 2007;97(6):1603–16. pmid:17304558
  35. 35. Mosahebi A, Woodward B, Wiberg M, Martin R, Terenghi G. Retroviral labeling of Schwann cells: In vitro characterization and in vivo transplantation to improve peripheral nerve regeneration. Glia. 2001;34(1):8–17. pmid:11284015
  36. 36. Guenard V, Kleitman N, Morrissey TK, Bunge RP, Aebischer P. Syngeneic Schwann cells derived from adult nerves seeded in semipermeable guidance channels enhance peripheral nerve regeneration. The Journal of Neuroscience. 1992;12(9):3310. pmid:1527582
  37. 37. Bunge RP. The role of the Schwann cell in trophic support and regeneration. Journal of Neurology. 1994;242(1):S19–S21. pmid:7699403
  38. 38. Colazzo F, Sarathchandra P, Smolenski RT, Chester AH, Tseng Y-T, Czernuszka JT, et al. Extracellular matrix production by adipose-derived stem cells: Implications for heart valve tissue engineering. Biomaterials. 2011;32(1):119–27. pmid:21074262
  39. 39. Hobson MI, Green CJ, Terenghi G. VEGF enhances intraneural angiogenesis and improves nerve regeneration after axotomy. Journal of Anatomy. 2000;197(4):591–605. pmid:11197533
  40. 40. Sanen K, Martens W, Georgiou M, Ameloot M, Lambrichts I, Phillips J. Engineered neural tissue with Schwann cell differentiated human dental pulp stem cells: potential for peripheral nerve repair? Journal of Tissue Engineering and Regenerative Medicine. 2017;11(12):3362–72. pmid:28052540
  41. 41. Ardakani AG, Cheema U, Brown RA, Shipley RJ. Quantifying the correlation between spatially defined oxygen gradients and cell fate in an engineered three-dimensional culture model. Journal of The Royal Society Interface. 2014;11(98):20140501. pmid:24966240
  42. 42. Fassnacht D, Pörtner R. Experimental and theoretical considerations on oxygen supply for animal cell growth in fixed-bed reactors. Journal of Biotechnology. 1999;72(3):169–84. pmid:10443023
  43. 43. Hay PD, Veitch AR, Smith MD, Cousins RB, Gaylor JDS. Oxygen Transfer in a Diffusion-Limited Hollow Fiber Bioartificial Liver. Artificial Organs. 2000;24(4):278–88. pmid:10816201
  44. 44. Croll TI, Gentz S, Mueller K, Davidson M, O’Connor AJ, Stevens GW, et al. Modelling oxygen diffusion and cell growth in a porous, vascularising scaffold for soft tissue engineering applications. Chemical Engineering Science. 2005;60(17):4924–34.
  45. 45. Popel AS. Theory of oxygen transport to tissue. Crit Rev Biomed Eng. 1989;17(3):257–321. pmid:2673661.
  46. 46. Dvorak HF. Vascular Permeability Factor/Vascular Endothelial Growth Factor: A Critical Cytokine in Tumor Angiogenesis and a Potential Target for Diagnosis and Therapy. Journal of Clinical Oncology. 2002;20(21):4368–80. pmid:12409337.
  47. 47. Carmeliet P, Jain RK. Molecular mechanisms and clinical applications of angiogenesis. Nature. 2011;473(7347):298–307. pmid:21593862
  48. 48. Aubert M, Chaplain MAJ, McDougall SR, Devlin A, Mitchell CA. A Continuum Mathematical Model of the Developing Murine Retinal Vasculature. Bulletin of Mathematical Biology. 2011;73(10):2430–51. pmid:21286832
  49. 49. Billy F, Ribba B, Saut O, Morre-Trouilhet H, Colin T, Bresch D, et al. A pharmacologically based multiscale mathematical model of angiogenesis and its use in investigating the efficacy of a new cancer treatment strategy. Journal of Theoretical Biology. 2009;260(4):545–62. pmid:19615383
  50. 50. Jain H, Jackson T. A Hybrid Model of the Role of VEGF Binding in Endothelial Cell Migration and Capillary Formation. Frontiers in Oncology. 2013;3(102). pmid:23675570
  51. 51. Maggelakis SA, Savakis AE. A mathematical model of growth factor induced capillary growth in the retina. Mathematical and Computer Modelling. 1996;24(7):33–41. https://doi.org/10.1016/0895-7177(96)00124-0.
  52. 52. Vempati P, Popel AS, Mac Gabhann F. Formation of VEGF isoform-specific spatial distributions governing angiogenesis: computational analysis. BMC Systems Biology. 2011;5(1):59. pmid:21535871
  53. 53. Wood BD, Quintard M, Whitaker S. Calculation of effective diffusivities for biofilms and tissues. Biotechnology and Bioengineering. 2002;77(5):495–516. pmid:11788949
  54. 54. Quintard M, Whitaker S. Coupled, nonlinear mass transfer and heterogeneous reaction in porous media. In: Vafai K, editor. Handbook of Porous Media. 2nd ed. Boca Raton, Florida, USA: Taylor & Francis Group, LLC; 2005. p. 3–37.
  55. 55. Tao J, Hu Y, Wang S, Zhang J, Liu X, Gou Z, et al. A 3D-engineered porous conduit for peripheral nerve repair. Scientific Reports. 2017;7(1):46038. pmid:28401914
  56. 56. Yin K, Divakar P, Hong J, Moodie KL, Rosen JM, Sundback CA, et al. Freeze-cast Porous Chitosan Conduit for Peripheral Nerve Repair. MRS Advances. 2018;3(30):1677–83. Epub 2018/02/20. pmid:30009044
  57. 57. Zhu L, Jia S, Liu T, Yan L, Huang D, Wang Z, et al. Aligned PCL Fiber Conduits Immobilized with Nerve Growth Factor Gradients Enhance and Direct Sciatic Nerve Regeneration. Advanced Functional Materials. 2020;30(39):2002610.
  58. 58. Weissberg HL. Effective Diffusion Coefficient in Porous Media. Journal of Applied Physics. 1963;34(9):2636–9.
  59. 59. Klemens PG. Thermal conductivity of composites. International Journal of Thermophysics. 1990;11(5):971–6.
  60. 60. Beeckman JW. Mathematical description of heterogeneous materials. Chemical Engineering Science. 1990;45(8):2603–10. https://doi.org/10.1016/0009-2509(90)80148-8.
  61. 61. Kokai LE, Lin Y-C, Oyster NM, Marra KG. Diffusion of soluble factors through degradable polymer nerve guides: Controlling manufacturing parameters. Acta Biomaterialia. 2009;5(7):2540–50. pmid:19369123
  62. 62. Ezra M, Bushman J, Shreiber D, Schachner M, Kohn J. Porous and Nonporous Nerve Conduits: The Effects of a Hydrogel Luminal Filler With and Without a Neurite-Promoting Moiety. Tissue Engineering Part A. 2016;22(9–10):818–26. pmid:27102571.
  63. 63. Meek MF, Den Dunnen WFA. Porosity of the wall of a Neurolac® nerve conduit hampers nerve regeneration. Microsurgery. 2009;29(6):473–8. pmid:19308952
  64. 64. Vleggeert-Lankamp CLAM, Wolfs J, Pêgo AP, Berg Rvd, Feirabend H, Lakke E. Effect of nerve graft porosity on the refractory period of regenerating nerve fibers. 2008;109(2):294. pmid:18671643
  65. 65. Oh SH, Kim JH, Song KS, Jeon BH, Yoon JH, Seo TB, et al. Peripheral nerve regeneration within an asymmetrically porous PLGA/Pluronic F127 nerve guide conduit. Biomaterials. 2008;29(11):1601–9. pmid:18155135
  66. 66. Aebischer P, Guenard V, Brace S. Peripheral nerve regeneration through blind-ended semipermeable guidance channels: effect of the molecular weight cutoff. The Journal of Neuroscience. 1989;9(10):3590. pmid:2795143
  67. 67. Cheema U, Rong Z, Kirresh O, MacRobert AJ, Vadgama P, Brown RA. Oxygen diffusion through collagen scaffolds at defined densities: implications for cell survival in tissue models. Journal of Tissue Engineering and Regenerative Medicine. 2012;6(1):77–84. pmid:21312340
  68. 68. Han P, Bartels DM. Temperature Dependence of Oxygen Diffusion in H2O and D2O. The Journal of Physical Chemistry. 1996;100(13):5597–602.
  69. 69. Gabhann FM, Ji JW, Popel AS. VEGF gradients, receptor activation, and sprout guidance in resting and exercising skeletal muscle. Journal of Applied Physiology. 2007;102(2):722–34. pmid:17038488.
  70. 70. Carreau A, Hafny-Rahbi BE, Matejuk A, Grillon C, Kieda C. Why is the partial oxygen pressure of human tissues a crucial parameter? Small molecules and hypoxia. Journal of Cellular and Molecular Medicine. 2011;15(6):1239–53. pmid:21251211
  71. 71. Kut C, Mac Gabhann F, Popel AS. Where is VEGF in the body? A meta-analysis of VEGF distribution in cancer. British Journal of Cancer. 2007;97(7):978–85. pmid:17912242
  72. 72. Robertson AM, Huxley C, King RHM, Thomas PK. Development of early postnatal peripheral nerve abnormalities in Trembler-J and PMP22 transgenic mice. Journal of Anatomy. 1999;195(3):331–9. pmid:10580849
  73. 73. Murray JD. Mathematical biology: I. An introduction: Springer Science & Business Media; 2007.
  74. 74. Gerhardt H, Golding M, Fruttiger M, Ruhrberg C, Lundkvist A, Abramsson A, et al. VEGF guides angiogenic sprouting utilizing endothelial tip cell filopodia. Journal of Cell Biology. 2003;161(6):1163–77. pmid:12810700
  75. 75. Jenq C-B, Coggeshall RE. Nerve regeneration through holey silicone tubes. Brain Research. 1985;361(1):233–41. pmid:4084796
  76. 76. Jenq C-B, Coggeshall RE. Permeable tubes increase the length of the gap that regenerating axons can span. Brain Research. 1987;408(1):239–42. pmid:3594212
  77. 77. Kim DH, Connolly SE, Zhao S, Beuerman RW, Voorhies RM, Kline DG. Comparison of Macropore, Semipermeable, and Nonpermeable Collagen Conduits in Nerve Repair. J Reconstr Microsurg. 1993;9(06):415–20. pmid:8283421
  78. 78. Huang L, Zhu L, Shi X, Xia B, Liu Z, Zhu S, et al. A compound scaffold with uniform longitudinally oriented guidance cues and a porous sheath promotes peripheral nerve regeneration in vivo. Acta Biomaterialia. 2018;68:223–36. pmid:29274478
  79. 79. Chamberlain LJ, Yannas IV, Arrizabalaga A, Hsu HP, Norregaard TV, Spector M. Early peripheral nerve healing in collagen and silicone tube implants: Myofibroblasts and the cellular response. Biomaterials. 1998;19(15):1393–403. pmid:9758039
  80. 80. Grützkau A, Krüger-Krasagakes S, Baumeister H, Schwarz C, Kögel H, Welker P, et al. Synthesis, Storage, and Release of Vascular Endothelial Growth Factor/Vascular Permeability Factor (VEGF/VPF) by Human Mast Cells: Implications for the Biological Significance of VEGF206. Molecular Biology of the Cell. 1998;9(4):875–84. pmid:9529385