Patient‐specific apparent diffusion maps used to model nutrient availability in degenerated intervertebral discs

Abstract

Introduction

In this study, magnetic resonance imaging data was used to (1) model IVD‐specific gradients of glucose, oxygen, lactate, and pH; and (2) investigate possible effects of covariate factors (i.e., disc geometry, and mean apparent diffusion coefficient values) on the IVD’s microenvironment. Mathematical modeling of the patient’s specific IVD microenvironment could be important when selecting patients for stem cell therapy due to the increased nutrient demand created by that treatment.

Materials and Methods

Disc geometry and water diffusion coefficients were extracted from MRIs of 37 patients using sagittal T1‐weighted images, T2‐weighted images, and ADC Maps. A 2‐D steady state finite element mathematical model was developed in COMSOL Multiphysics® 5.4 to compute concentration maps of glucose, oxygen, lactate and pH.

Results

Concentration of nutrients (i.e., glucose, and oxygen) dropped with increasing distance from the cartilaginous endplates (CEP), whereas acidity levels increased. Most discs experienced poor nutrient levels along with high acidity values in the inner annulus fibrosus (AF). The disc’s physiological microenvironment became more deficient as degeneration progressed. For example, minimum glucose concentration in grade 4 dropped by 31.1% compared to grade 3 (p < 0.0001). The model further suggested a strong effect of the following parameters: disc size, AF and CEP diffusivities, metabolic reactions, and cell density on solute concentrations in the disc (p < 0.05).

Conclusion

The significance of this work implies that the individual morphology and physiological conditions of each disc, even among discs of the same Pfirrmann grade, should be evaluated when modeling IVD solute concentrations.

INTRODUCTION

Back pain is the leading cause of disability worldwide, causing a significant socioeconomic impact.¹ A source of back pain can be the degeneration of the intervertebral disc (IVD), allowing nerve ingrowth, facet joint arthritis, and disc bulging or osteophyte formations that press on nearby nerve roots or the spinal cord.^2,3 Animal and human studies have shown promising potential for cell injections to restore IVD health. In a pilot study performed on 14 patients, chondrocyte injections demonstrated the ability to maintain disc height and hydration level.^4,5 The results showed a significant reduction in the disability index in treated patients 2 years post discectomy.^4,5 Introducing disc cells into the annulus fibrosus (AF) and nucleus pulposus (NP) of animal models maintained the demarcation between the AF and the NP. They also delayed the degeneration process of the IVD by increasing the remodeling of the extracellular matrix (ECM).^6-8 Studies assessing the viability of transplanted cells showed good cellular viability and proliferation rates after treatment.^8,9 In rabbit models, transplanted human NP cells increased disc height, the expression of aggrecan, and type II collagen.¹⁰ The results of these studies implied good integration of injected cells into the disc's environment, which demonstrated the promising potential in using cell therapy to regenerate the disc.^4,6-11

There are, however, limitations to relying on cell therapy to restore the IVD due to the morphological and biochemical changes attributed to IVD degeneration. First, degenerated IVDs suffer from a loss in their proteoglycan content, leading to IVD dehydration. Bulk water content is critical for IVD health as it contributes to molecular transport and the IVD's biomechanics.^12,13 Second, structural changes of the AF and the cartilaginous endplates (CEP) limit molecular transport through the IVD.^12,13 Lastly, degeneration of the IVD leads to changes in the NP's ECM that altered diffusion kinetics of solutes through the IVD.^14-16 As a result of the changes in IVD biology, diffusion rates decrease dramatically, leading to poor nutrition (ie, glucose and oxygen), high acidity, and high cell death.^17-30

In addition, clinical stem cell studies reported widely varying cell doses being administered to treat IVD degeneration, ranging from 1 million to 40 million cells per patient as reviewed by Sakai and Schol.^{4,21,22,31-38} Considering the lack of nutrient supply in the human IVD,¹³ we believed it was essential to model the patient-specific nutrient distributions before cell therapy and identify the appropriate dose of cell injections based on the nutritional capacity of the individual disc.³⁷

Mathematical models used experimentally measured transport properties, including diffusion coefficients for glucose, oxygen, and lactate.^14,26,39-43 In addition to IVD properties such as IVD water content, cell density, and the rate of cellular glycolysis as a function of oxygen content to simulate solute transport phenomena in the disc.^14,26,39-43

These models predicted a high concentration of nutrients in the outer regions of the disc: closer to the endplates and the outer AF, compared to the mid-axial center of the disc and in the inner AF.^44-46 The results from said models agreed with experimentally measured solute concentrations in the disc. Oxygen concentration profiles measured during discography were high in the outer regions of the disc and gradually decreased toward the mid-axial center of the disc.⁴⁷ Lactate concentrations measured in the NP tissue taken by biopsy from patients and excised segments of the disc were low in the outer regions of the disc and higher in the mid-axial center of the disc.^47,48 There was no literature data about actual glucose measurements in human IVDs. Since glucose is required to produce lactate, their concentrations are necessarily coupled.⁴⁶ Therefore, glucose levels should decrease toward the mid-axial center of the disc.⁴⁹ Lastly, studies measured pH values in the NP of patients with sciatica to range between 5.7 and 7.5.^48,50

This study aimed to use data from magnetic resonance imaging (MRI) to model patient-specific IVD distributions of glucose, oxygen, and lactate and investigate factors affecting solute availability in the disc. Disc geometry was extracted from MR T1-weighted (T_1w) and T2-weighted images (T_2w). NP diffusion coefficients were calculated from apparent diffusion coefficient (ADC) maps. The first part of this study, looked at significant differences in solute concentrations between the different degeneration grades. We hypothesized the following with increasing IVD degeneration: (1) average nutrient concentration in the entire disc would decrease, accompanied by an increase in lactate concentrations; (2) minimum concentration values of nutrients would decrease accompanied by an increase in maximum lactate concentrations (3) average solute concentrations in the disc would correlate with their extreme solute values (minimum/maximum) in the inner AF. Then we looked at the effect of disc area, disc height, and average ADC on solute concentrations.

METHODS Patient MR data

This study was approved by the Institutional Ethics Committee at the Irkutsk Scientific Center of Surgery and Traumatology, Irkutsk, Russia (identifier no. 15-15-30 037). All patients voluntarily provided their written consent to take part in this study. A total of 37 IVDs were randomly chosen from a preexisting MR patient database.²⁷ Each Pfirrmann grade group (grades 2-4) had 10 randomly selected MR scans of IVDs. Only seven sets of IVD MR scans with degeneration grade 5 were available for modeling; therefore, randomization was not possible for this group. The IVDs included in this study represented disc levels between L3L4 and L5S1 due to the wrap-around artifacts that distorted higher levels, that is, L1L2 and L2L3. No MR scans of grade 1 were available for modeling since they suffered from image artifacts. Patient characteristics were summarized in Table 1.

TABLE 1 Summary of patient characteristics

Patient characteristic	Grade 2 (n = 10)	Grade 3 (n = 10)	Grade 4 (n = 10)	Grade 5 (n = 7)	P-valuê
Age	36.5 ± 7.0	36.6 ± 6.6	37.4 ± 7.5	48.4 ± 6.8	0.0014
BMI	142.2 ± 51	101.5 ± 37	92.8 ± 17	103 ± 29	0.4770
Height (cm)	116.8 ± 49	152 ± 47	175 ± 6.4	163 ± 45	0.2534
Sex (m:f)	9:1	9:1	10:0	7:0	0.6442
Disc pathology	—	—	—	—	1.82E−06
None	90.0%	80.0%	20.0%	0.0%	—
Herniation	0.0%	20.0%	70.0%	57.1%	—
Spondylosis	0.0%	0.0%	10.0%	28.6%	—
Stenosis	10.0%	0.0%	0.0%	14.3%
Disc level	—	—	—	—	0.0241
L3L4	70.0%	30.0%	10.0%	0.0%	—
L4L5	10.0%	60.0%	40.0%	28.6%	—
L5S1	20.0%	10.0%	50.0%	71.4%	—

Note: Ten IVDs from each degeneration grade were randomly chosen from the 21 patients. Randomization could not be applied for grade 5 because only 7 IVDs were available. IVDs included in this study represent levels between L3L4 and L5S1 due to image artifacts that distorted higher levels, that is, L1L2 and L2L3. The table also consists of each patient's age and IVD pathology condition. ̂P-values reported are from one-way ANOVA or Kruskal-Wallis when data were parametric or nonparametric, respectively.

Magnetic resonance imaging

Briefly, Sagittal T1-weighted (T_1w), T2-weighted (T_2w), and diffusion weighted-imaging (DWI) scans of human lumbar discs were obtained using a 1.5 T Siemens Magnetom Essenza scanner (Siemens Healthineers, Erlangen, Germany).²⁷ Images had a 30 × 30 cm field of view. The following matrix sizes 320 × 320, 384 × 384 were used for T_1w, and T_2w, with pixel sizes 1.0438 × 1.0438 (mm/pixel) and 0.8698 × 0.8698 (mm/pixel), respectively.²⁷ DWIs were collected at three different b-values (b = 50, 400, 800 seconds/mm²) using a body coil with a TR of 3000 ms, TE of 93 ms, six averages, a matrix size of 156 × 192, and a pixel size of 1.7396 × 1.7396 (mm/pixel).²⁷ The mid-sagittal slice of each MR imaging modality was used to generate 2D models of the disc. Time of day when the scans were performed was not controlled for.²⁷

Pfirrmann grading

Discs were evaluated for degeneration and assigned a Pfirrmann grade by three neurosurgeons specializing in spinal surgery with 5 to 15 years of experience.^27,51 Each surgeon graded all discs then collectively discussed and formed a consensus on any conflicting grades.

Segmentation of disc masks

Mask segmentation was performed in Materialize Mimics 21.0 (Leuven, Belgium) (Figure 1). T_1w images were used to segment CEP masks by hand. Whole IVD and NP masks were segmented using histogram-based thresholding on T_2w images. AF masks were generated by subtracting NP and CEP masks from complete IVD masks. The masks were cropped in MATLAB 9.5 (Natick, Massachusetts) to show the disc of interest. Further processing of masks, which included manual subtraction of masks from one another, was performed to prevent overlap between masks (in most cases, this last step was not needed).

FIGURE 1. Annulus fibrosus (AF), cartilaginous endplates (CEP), and nucleus pulposus (NP) masks (in white left to right) were segmented in Materialize Mimics 21.0 for each patient. AF, CEP, and NP masks from intervertebral disc #1 L3L4 were segmented then overlaid over a T2w scan in this figure

Generation of ADC map

Series of DWI scans were processed in MATLAB 9.5 to generate ADC maps for the NP using the mid-sagittal slice (Figure 2).

FIGURE 2. Diffusion weighted-imaging (left) for intervertebral disc #1, level L3L4, unprocessed apparent diffusion coefficient (ADC) map (middle), cropped ADC map (right)

ADC values were calculated using the Stejskal-Tanner equation as established by Le Bihan (Equation 1).⁵² [Image Omitted. See PDF]Where Equation (1), $S_0$ was the MR signal intensity at baseline, $S$ the MR signal intensity after diffusion gradients had been applied, and $b$ was the attenuation factor that depended on the gradient pulse. Raw ADC maps were further processed by replacing both negative and NAN (not a number) values with “0.” Finally, ADC maps were cropped to match the dimensions of the corresponding NP mask. The cropped maps were used in the model to compute diffusion in the disc.

Model transport properties

A two-dimensional steady-state mathematical model was developed in COMSOL Multiphysics 5.4 (Stockholm, Sweden). The temperature was set to $37°\mathrm{C}$ and diffusion throughout the disc was assumed to be isotropic. The simulated time in the model was 30 hours, and it converged with relative tolerance of 0.5%.

Model geometry

AF, CEP, and NP masks were overlaid on top of a square (45 mm × 45 mm) which represented a well-mixed container in which solute transport took place. Disc masks were scaled to their true size and always positioned at the center of the geometric shape. A physics-controlled mesh with extra fine element size (maximum element size 0.9 mm, minimum element size 0.00338 mm, curvature factor 0.25, maximum element growth rate 1.2) was implemented with Adaptive Mesh Refinement (Figure 3).

FIGURE 3. Adaptive mesh refinement for intervertebral disc #1, level L3L4. COMSOL Multiphysics had automatically inserted additional mesh elements to minimize error, increase accuracy. and decrease computational time (left). Zoomed view of mesh boundary layers (right) shows the areas around the IVD boundaries that contained additional smaller mesh elements

Boundary and initial conditions

Boundary concentration values for glucose, oxygen, and lactate (5 (mol/m³), 0.06 (kPa), 0.9 (mol/m³), respectively) were chosen from previously published IVD finite element models.^42,44,45,47 Boundary conditions were assigned to the sides of the geometric shape to simulate a disc suspended in a well-mixed solution. The assumption of a well-mixed solution was essential to eliminate variability in solute gradients outside the disc between different models (Figure 4). Physiological levels of glucose, oxygen, and lactate were assumed in the well-mixed solution.^44-46 Initial conditions inside the disc for glucose, oxygen, and lactate were set to 1 (mol/m³), 0.1 (kPa), and 0 (mol/m³), respectively.^44-46

FIGURE 4. Glucose diffusion map of intervertebral disc #1, disc level L3L4. Incorporating apparent diffusion coefficient (ADC) maps to calculate diffusion in the nucleus pulposus reflected the local variation in diffusion coefficients. No gradient could be observed in the cartilaginous endplates and annulus fibrosus because ADC maps were not obtained at a resolution high enough to allow for extraction of diffusion coefficients

Constant parameters

Parameter maps for solute diffusion, water content, and cell density for each segment of the disc were generated in COMSOL Multiphysics by multiplying the mask with the corresponding parameter value⁴² (Figure 4). Discs with degeneration grade of 2 to 3 were considered healthy, whereas 4 to 5 grade discs were considered degenerate (Table 2).

TABLE 2 The model's constant parameters

IVD region	Water content (%)		Diffusion coefficient (m2/s)			Cell density (1/m3)
—	Healthy disc	Degenerate disc	Glucose	Oxygen	Lactate	—
NP	83[53]	78[53]	see “Section 2.10”			4.0E+12[49]
AFa	69[54]	53[54]	2.85E−10[42]	1.1E−9[42]	4.24E−10[42]	9.0E+12[49]
CEP	60[45]	42[45]	2.11E−10[42]	7.81E−10[42]	3.14E−10[42]	1.5E+13[49]

Note: IVD parameters for water content, diffusion coefficients for glucose, oxygen, and lactate, and cell density were taken from other published models. These parameters were multiplied with each mask to create an IVD map of varying parameter values based on each anatomical part.

^aAverage values of inner and outer AF.

Diffusion coefficients in NP

Solute diffusion coefficients in the NP were generated by converting water diffusion coefficients from ADC maps to the corresponding solute diffusion coefficient. Factors used to scale diffusion coefficients of water to glucose, oxygen, and lactate were 0.319, 1.02, and 0.478, respectively. These factors were generated by scaling the value of water diffusion in tissue to solute diffusion in tissue using hindered solute diffusion theory as developed by Renkin.⁵⁵ His equation describes the diffusion of a solute molecule (ie, glucose, oxygen, and lactate) through a tiny capillary pore filled with a liquid solvent (ie, water). As the size of the solute increased, the diffusive transport through the solvent was hindered by the presence of the pore, specifically the pore wall.^55,56 Renkin's model was used to describe diffusion in the IVD because the ECM's pore diameter is 7 to 20 times the molecular diameter of glucose, oxygen, and lactate.⁵⁷

Briefly, diffusion of molecules through porous material could be characterized using two controlling mechanisms, geometric tortuosity and hindered diffusion. Geometric tortuosity represents the ratio between the average pathways and the straightest one taken by a specific molecule as it traverses from point A to B in a porous medium.^58-60 Hindered diffusion describes the molecular interactions with the pore's walls, which inhibit solute transport as described by Renkin's equation. In situations where the pore size is over 30 to 40 times the diameter of molecules, diffusion is controlled by the geometric tortuosity theory.⁶⁰ When the porous medium's pore size is 10 times or less the molecule's diameter, then hindered diffusion is the dominating factor.⁶⁰ This hindrance was modeled using Equation (2).⁵⁵ [Image Omitted. See PDF]In this equation both correction factors, $F_1$ and $F_2$, are functions of the reduced pore diameter $\left(\phi \right)$, and are theoretically bounded by 0 and 1 as shown in Equation (3).⁵⁵ [Image Omitted. See PDF] $F_1\left(\phi \right)$ is known as the stearic partition coefficient and is based on geometric arguments for stearic exclusion in Equation (4).⁵⁵ [Image Omitted. See PDF]The correction factor $F_2\left(\phi \right)$ is known as the hydrodynamic hindrance factor. It is based on several hydrodynamic calculations, including the hindered Brownian motion of the solute within the solvent-filled pore.⁵⁵ Renkin developed the following relationship for $F_2\left(\phi \right)$, assuming the solute is a rigid sphere diffusing through a straight cylindrical pore (Equation 5).⁵⁵ [Image Omitted. See PDF]

Based on this analysis, ADC values were scaled for a given solute based on a ratio of literature values for solute diffusion in water and the diffusion of water in water. It was assumed that all four solutes (glucose, oxygen, lactate, and water) were sufficiently small (on the scale of 300-900 pm) relative to the pore to approximate $F_1\left(\phi \right)$ and $F_2\left(\phi \right)$ as the same across all solutes.

The diffusion of water in tissue was modeled as hindered-self diffusion of water through a solvent-filled pore, in which diffusion constant parameters $F_1\left(\phi \right)$ and $F_2\left(\phi \right)$ were already accounted for in the ADC values. Based on this method, scaling factors were generated to relate water diffusion in the tissue (extracted from ADC maps) to solute diffusion of glucose, oxygen, and lactate in the tissue.

A comparison of the diffusion coefficients used in this model and the reported values in the literature showed that this model's diffusion coefficients fell within the range of the reported values (Table 3).

TABLE 3 Comparison of this model's diffusion coefficients and literature values

Year	Author	Model	IVD region	Dg (m2/s)	Dl (m2/s)	Do (m2/s)
This model			CEP	6.73E−11	1.50E−10	7.96E−10
			AF	9.08E−11	2.02E−10	1.07E−09
			NP	1.22E−10 to 2.54E−10	1.83E−10 to 3.81E−10	3.90E−10 to 8.13E−10
2016[30]	Cisewski	FE model	CEP	2.68E−11	4.52E−11	—
2016[30	Wu	Human ex vivo	CEP	2.68E−11	4.52E−11	—
2012[61]	Jackson	Human in vitro	AF	3.56E−11 to 8.71E−11	—	1.13E−10 to 1.85E−9
2009[62]	Jackson	Bovine in vitro	AF	—	—	1.56E−09
2008[63]	Jackson	Bovine in vitro	AF	1.38E−10 to 9.17E−11	—	—
2009[64]	Magnier	FE model	CEP	9.17E−10	1.39E−09	3.00E−09
2007[42]	Soukane	FE model	CEP	2.11E−10	3.14E−10	7.81E−10
			AF	2.85E−10	4.24E−10	1.05E−09
			NP	3.78E−10	5.61E−10	1.39E−09
1975[49]	Maroudas	Human in vitro	CEP	2.43E−10	—	—
			AF	2.50E−10	—	—
2003[65]	Sélard	FE model	AF	2.50E−10	4.86E−10	8.33E−10
			NP	3.75E−10	6.11E−10	1.28E−09

Note: Diffusion coefficients of glucose, oxygen, and lactate used in this model were compared with their corresponding values reported in the literature. This model's coefficients, which were extracted from ADC maps and multiplied by scaling factors, fell within the range of diffusion coefficients reported in the literature. The range of diffusion coefficients used in the NP for this model depicted the highest and lowest values obtained from the ADC maps.

Governing equations

Diffusion of solutes (glucose, oxygen, and lactate) was modeled by substituting Fick's first law for flux into the equation of mass conservation where $C$ was solute concentration (mol/m³), $t$ was time (s), $D$ was solute diffusion coefficient (m²/s), $\dot{R}$ was the reaction term for either consumption or generation of solutes (mol/m³ s), and $x$ was diffusion distance (m) (Equation 6).⁶⁶ [Image Omitted. See PDF]The reaction term in Equation (6) was represented by nonlinear coupled equations for consumption of oxygen (nmol/million cells h) (Equation 7) and generation of lactate (nmol/million cells h) (Equation 8).³⁹ [Image Omitted. See PDF] [Image Omitted. See PDF]Oxygen solubility in water was used to convert oxygen consumption levels to kPa via conversion factors (Equation 9).⁴⁶ [Image Omitted. See PDF]pH values were correlated to lactate concentrations in the disc adapted from a study that investigated the correlation between lactate levels and pH in patient discs with lumbar sciatica (Equation 10).⁴⁸ [Image Omitted. See PDF]Glucose consumption was predicted based on lactate production. Bibby et al estimated the ratio of lactic production to glucose consumption to 2.01, which was expected from glycolytic pathways (Equation 11).³⁹ [Image Omitted. See PDF]Reaction equations were multiplied by water content and cell density to arrive at a final reaction unit of (mol/m³ s).

Calculation of total solute mass and average concentration

Average operation function in COMSOL multiphysics (Equation 12) was used to calculate average solute concentration (mol/m²) within the boundaries of the disc.[Image Omitted. See PDF]where, aveop# was the average operation function, and IVD(x,y) defined the boundaries of the disc.

Sensitivity analysis

A sensitivity analysis (SA) was performed on three randomly selected models of healthy IVDs to find the parameters with the most significant effect on nutrient distribution. The following model variables were changed by ±50% from baseline values: diffusivity of the AF, CEP, and NP, reaction rates for glucose, lactate, oxygen, disc size, and NP cell density. This SA was done to determine parameters that significantly affected both average concentrations in the IVD and minimum and maximum concentrations in the inner AF. A general regression model was performed on the % change in solute concentrations to find the effect size of model variables (P <.05).

Statistical analysis

For this study, 10 discs were randomly selected to represent Pfirrmann degeneration grades 2 through 5 from a pre-existing patient MR database. Randomization could not be applied for degeneration grade 5 because only 7 IVD datasets were available. Discs included in this study represented levels between L3L4 and L5S1 due to the wrap-around artifacts that distorted higher levels, that is, L1L2 and L2L3 (Table 1). Sample size in this study provided a statistical power of 95% for the mean difference of the following modeled parameters between different Pfirrmann grades: minimum glucose concentration, minimum oxygen concentration, maximum lactate concentration, average glucose concentration and average oxygen concentration. Average lactate concentration had a statistical power of 42% due to significant variance.

The homogeneity assumption of variance was verified using Levene's test (P <.05). Analysis of covariance (ANCOVA) was performed in TIBCO Statistica v13.5 (Palo Alto, California) to find a significant difference in mean solute concentrations and mean minimum and maximum solute concentrations between disc degeneration grades. ANCOVA was followed with a post-hoc Tukey honestly significant difference (HSD) test. Disc height, disc area, and ADC mean intensity were covariate factors. Outliers were measured using inclusive median quartile calculation. Partial correlations were performed between the estimated metrics (average solute concentrations in the disc as well as minimum and maximum solute concentrations in the inner AF) and disc properties (disc area, disc height, and average ADC intensity). Statistical analysis was performed on patient demographics using one-way ANOVA or Kruskal-Wallis (P <.05) when data were parametric or nonparametric, respectively. Finally, one-way ANOVA followed by a post-hoc test was performed to investigate the possible effects of patient characteristics, that is, sex, IVD pathology, and IVD level on solute concentrations.

RESULTS Nutrient distribution

For all degeneration grades, the concentration of nutrients (ie, glucose and oxygen) dropped with distance from the CEP. In contrast, acidity levels increased especially toward the center of the disc (Table 4). Most discs experienced the lowest nutrient levels along with the highest lactate in the inner AF (Table 4). The model predicted minimum glucose concentrations to range between 0.393 (mol/m³) in most degenerated IVDs to 3.43 (mol/m³) in least degenerated IVDs. Minimum oxygen and pH values ranged from 4.67 (kPa) and 7.18 in least degenerated IVDs to 2.76 (kPa) and 6.96 in most degenerated IVDs.

TABLE 4 Examples of solute and pH distribution models of discs with different degeneration grades

Patient	Grade level	T2w	ADC	Glucose gradient (mM)[IMAGE OMITTED. SEE PDF.]	Lactate gradient (mM)[IMAGE OMITTED. SEE PDF.]	Oxygen gradient (kPa)[IMAGE OMITTED. SEE PDF.]
1	2	[IMAGE OMITTED. SEE PDF.]	[IMAGE OMITTED. SEE PDF.]	[IMAGE OMITTED. SEE PDF.]	[IMAGE OMITTED. SEE PDF.]	[IMAGE OMITTED. SEE PDF.]
4	3	[IMAGE OMITTED. SEE PDF.]	[IMAGE OMITTED. SEE PDF.]	[IMAGE OMITTED. SEE PDF.]	[IMAGE OMITTED. SEE PDF.]	[IMAGE OMITTED. SEE PDF.]
19	4	[IMAGE OMITTED. SEE PDF.]	[IMAGE OMITTED. SEE PDF.]	[IMAGE OMITTED. SEE PDF.]	[IMAGE OMITTED. SEE PDF.]	[IMAGE OMITTED. SEE PDF.]
10	5	[IMAGE OMITTED. SEE PDF.]	[IMAGE OMITTED. SEE PDF.]	[IMAGE OMITTED. SEE PDF.]	[IMAGE OMITTED. SEE PDF.]	[IMAGE OMITTED. SEE PDF.]

Note: All discs were presented with the posterior aspect on the left and the anterior side on the right. Each model was constructed in COMSOL Multiphysics 5.4 using the middle sagittal slice of each patient's T_2w, T_1w, and DWI scans. pH distributions (linearly correlated to lactate levels) were not shown. The lowest levels of nutrients and the highest acidity levels appeared at the inner AF region in most discs (refer to Table S1 for complete model results).

Between Pfirrmann grades

There was a significant difference in patient demographics, including age (P = .0014), IVD level (P = .0241), and IVD pathology (P = 1.82E−06), between the different grades (Table 1). The model demonstrated that the IVD's nutrient environment became harsh as degeneration progressed. For example, average glucose concentrations followed a decreasing trend that was more obvious in minimum glucose levels (Figure 5A). Oxygen levels also decreased as the disc degenerated and showed a similar trend to glucose (Figure 5C). The decrease in nutrient levels was accompanied by increased maximum lactate and acidity levels (Figure 5B2).

FIGURE 5. Estimation of concentration averages with SD in intervertebral discs (IVDs) of Pfirrmann grades 2–5. Average solute concentrations were calculated using COMSOL Multiphysics for (A1) glucose (mM), (B1) lactate (mM), (C1) oxygen (kPa). COMSOL Multiphysics identified critical solute values in the inner AF for (A2) minimum glucose (B2), maximum lactate (C2) minimum oxygen (* = 0.05, § = 0.005, # = 0.0001). Patients #5, 7, and 18 were outliers. Patients #14 and 20 have larger and small IVD areas, respectively (Figure 7)

Quantitatively, the model showed a more substantial effect of degeneration on the minimum and maximum solutes compared to average solute concentrations in the IVD. For example, as the disc degenerated from grade 3 to 4, the average lactate concentration dropped by 7.59% (Table 5). In contrast, maximum lactate levels increased by 28.2%. In addition, minimum glucose and oxygen levels decreased with degeneration by 31.1% and 13.9% between grades 3 and 4; whereas, average glucose and oxygen decreased by 2.54% and 1.53%, respectively (Table 5).

TABLE 5 Quantified changes in solute concentrations between different degeneration grades

Degeneration progression	Glucose	Lactate	Oxygen
Grade 2—Grade 3	−2.79%	4.84%	−0.91%
Grade 3—Grade 4	−2.54%	7.59%	−1.53%
Grade 4—Grade 5	1.55%	−3.81%	1.10%
Degeneration progression	Minimum glucose	Maximum lactate	Minimum oxygen
Grade 2—Grade 3	−4.80%	4.79%	−2.17%
Grade 3—Grade 4	−31.1%	28.2%	−13.9%
Grade 4—Grade 5	16.5%	−8.00%	6.39%

Note: COMSOL Multiphysics 5.4 automatically identified minimum and maximum solute values in the IVD. However, all points fell within the inner AF (except for the following patients: #13, #15, and #18 in grade 4; and patient #20 in grade 5). (n = 10 for grades 2, 3, 4 and n = 7 for grade 5).

Between patients

The effects of patient demographics, that is, IVD level, IVD condition, and patient sex on solute concentrations were analyzed (Figure 6). Significant difference in average nutrient levels between IVD levels L3L4 and L4L5 (ANOVA, P <.05) (Figure 6A1). There was also a significant difference between herniated and none pathologic IVDs (ANOVA, P <.05) (Figure 6B2). Patient sex did not affect the mean solute concentrations, although only two IVDs from the same female patients were used in this study (Figure 6C). Mean patient age of grade 5 was 13 years higher than the mean patient age in all other groups (ANOVA, P <.005; Table 1). There was no strong evidence to support the variation in mean patient BMI and mean patient height between degeneration grades (ANOVA, P = .4770, and P = .2534, respectively; Table 1). No evidence was found to support a difference between patient sex between degeneration grades (Kruskal-Wallis, P = .6442; Table 1). Lastly, grade 2 has significantly more L3L4 IVDs compared to all other degeneration groups (Kruskal-Wallis, P <.05; Table 1).

FIGURE 6. Estimation of solute concentration means for glucose (mM), oxygen (kPa), and lactate (mM) in the IVD vs patient characteristics. (A1, A2) disc level, (B1, B2) IVD pathology condition, and (C1, C2) patient sex. An * indicated significance with a P-value of .05

The model further suggested a strong effect (ANCOVA, P <.05) between solute concentrations in the disc and the following three covariate factors (which were made patient-specific), disc height, disc area, and mean NP ADC values (Table 6). Average concentrations of all solutes correlated with all minimum and maximum values (P <.05; Table 6). Nutrient distributions varied significantly when visually compared model results in grades 3 to 5 (Tables 4 and S1).

TABLE 6 Partial correlations for different measured metrics in the disc

Variable	Min. glucose	Max. lactate	Min. oxygen	Ave. glucose	Ave. lactate	Ave. oxygen	Area	Height	Mean ADC intensity
Min. glucose	1.000	−1.000	0.997	0.789	−0.868	0.866	−0.602	−0.400	0.562
Max. lactate	−1.000	1.000	−0.997	−0.788	0.866	−0.864	0.604	0.407	−0.574
Min. oxygen	0.997	−0.997	1.000	0.769	−0.848	0.849	−0.608	−0.413	0.574
Ave. glucose	0.789	−0.788	0.769	1.000	−0.976	0.972	−0.670	−0.611	0.298
Ave. lactate	−0.868	0.866	−0.848	−0.976	1.000	−0.996	0.709	0.544	0.495
Ave. oxygen	0.866	−0.864	0.849	0.972	−0.996	1.000	−0.743	−0.607	−0.469

Note: Reported numbers are Pearson R correlation values of fitted linear regression between groups. The sign indicates a positive or negative correlation. All correlations were considered significant with P <.05. Grey highlighting indicates p<0.05.

Examining the outliers in this study showed that IVDs with a large IVD area and a high Pfirrmann score (grades 4 or 5) were more likely to have very extreme levels of nutrients and lactate, as was the case for patient #18 (Figures 5 and 7, Table S1). Patient #20 in grade 5 had a minimal IVD area due to a herniation, but solute levels fell within the expected values for that group (Figure 7, Table S1). Large healthy discs did not seem to suffer from nutrient deprivation. For instance, patient #14 in grade 3 had a large disc area (33% larger than the mean area for grade 3 discs), but the IVD's nutrient levels were within one SD from that group's mean (Figures 5 and 7, Table S1).

FIGURE 7. Average intervertebral disc area in pixels for each degeneration group. Outliers are marked with the corresponding patient number. An * indicated significance with a P-value of .05. Patients #14, 18, and 20 were outliers

The SA for this model showed that modifications to the IVD size had the most notable effect on all solute concentrations (Figure 8). The results showed an increase in IVD size to correlate with poor nutrient availability, while a decrease in size led to an increase in nutrient levels. For example, increasing the IVD size by 50% decreased minimum glucose by 80%. In contrast, a 50% decrease in IVD size resulted in a 63% increase in glucose concentration. A decrease in AF and CEP diffusivities lowered minimum glucose levels in the disc by 43% and 17%, respectively, while increased lactate concentrations by 36% and 9%. The NP's diffusivity affected maximum lactate levels only (Figure 8B). An increase in the consumption rates of glucose and oxygen showed a +10% drop in their minimum concentrations (Figure 8A,C) while it did not affect lactate levels (Figure 8B). Lactate production rates correlated with a +35% change in minimum glucose, and +30% change in maximum lactate concentrations (Figure 8A,B). More factors affected minimum glucose by +40% compared to other solutes, including disc size, glucose consumption rate, lactate production rate, AF, and CEP diffusivities (Figure 8A). Changes to the IVD's cell density affected maximum lactate and oxygen concentrations (Figure 8B,C).

FIGURE 8. Sensitivity analysis on model variables including, diffusivity of the annulus fibrosus, cartilaginous endplates, and nucleus pulposus, reaction rates (RXN) of glucose, lactate, and oxygen, cell density, and disc size. Solute % change from baseline values are presented for (A) glucose, (B) lactate, and (C) oxygen. General linear model was performed to find significance indicated by * (P = .05)

DISCUSSION

The main objective of this study was to create finite element models that incorporated patient-specific MRI data to model nutrient and metabolite distributions in IVDs of different Pfirrmann grades. The study also assessed the effects of disc area, disc height, and mean ADC on the IVD's nutrition (Table 6). The model estimated average glucose, oxygen, and lactate distributions in the disc, and minimum glucose, minimum oxygen, and maximum lactate concentrations. Taken together, these metrics identify potential low-nutrient regions in each disc, where cell growth is likely to be inferior.

Nutrient distribution

The model showed nutrient levels moving from physiological levels near the IVD's boundary to their minimums near the inner AF (Tables 4 and S1). Critical solute levels were expected in this part of the IVD due to the long diffusion pathway that solutes need to take from the CEP (Figure 9). Multiple studies have shown the dependency of nutrient availability on diffusion through the CEP, attributed to its porosity and low thickness.^29,67 The outer AF relies on molecular transport through its well-aligned collagen rings for nutrient supply.⁶⁷ Nonetheless, due to the AF's anisotropic properties and thicker layers, molecular transport was still less effective than through the CEP.⁶⁷ These IVD traits create regions that are vulnerable to nutrient deprivation, and are believed to not support the cellular nutrient demand (Table 4). The amount of glucose that is consumed in each degeneration grade was equal to about half of the amount of lactate being produced, when comparing these values to the boundary conditions. For example, in grade 2, the amount of consumed glucose was 1 (mM) on average (Figure 5A1) and the corresponding amount of lactate was about 1.8 (mM) on average (Figure 5B1). When comparing the amount of lactate being produced in degenerated grades, one can notice that it is higher for grades 4 and 5 (about 2 mM) compared with healthy IVDs. This increase cannot be attributed to changes in metabolic rates or cell density because both variables were held constant. Instead, it is attributed to the disrupted balance between metabolic rates and diffusion rates of lactate. The decrease in the ADC values and diffusion coefficients for high Pfirrmann grade IVDs, reduces the disc's ability to clear lactate at an efficient rate. In-vivo measurements and finite element models have already established these solute distribution characteristics in the disc. For example, oxygen levels in discs from patients with scoliosis were measured to be lower than plasma levels and progressively decrease toward the axial center of the disc, with minimums in the inner AF ranging between 0.67 (kPa) and 20 (kPa),⁴⁷ which matches our model's predictions (Tables 4 and S1, Figure 5C). In addition, several finite element models have been used to illustrate glucose distribution in human IVDs where they predict glucose levels to decrease from plasma levels at the outer regions of the disc to ~1 (mol/m³) in the inner regions of the disc. As shown in Figure 5A and Table S1, the model's estimation for glucose levels reflected literature values and followed similar distribution patterns to that of oxygen (Table S1).^{42,45,46,64,68,69} This model's prediction also agreed with measured lactate levels in human discs, which increased from 1 (mol/m³) (plasma level) to 6 (mol/m³) at the center of the disc (Figure 5B).⁴⁷ In this model, lactate was linearly correlated to pH, so regions with high lactate correspond to low pH values (Equation 10). Modeled pH levels fell within literature values of 6.7 to 7.4.^44-46 Furthermore, average glucose, oxygen, and lactate concentrations significantly correlated with their critical levels in the disc (Table 6).

FIGURE 9. Streamlines show lactate flux through the different parts of the IVD (annulus fibrosus [AF]: light gray, cartilaginous endplates [CEP]: dark gray, nucleus pulposus [NP]: white). The color scale reflects flux magnitude of lactate. High flux was observed near both CEPs, while the AF experienced low flux indicating better transport through the former. Long diffusion pathways extending from the inner AF through the CEP lead to low nutrient levels in the inner AF

Solute flux estimations reflected the critical role of the CEP in solute transport. Lactate flux density was greater in the CEP than the AF, indicating that the CEP was the path of least resistance in terms of molecular transport (Figure 9). In addition, predictions showed that solutes need to travel longer distances to reach or clear the inner AF (Figure 9). This circumstance may prevent nutrients from being replenished at an efficient rate; while lactate tends to accumulate and lead to high acidity.^12,25,70-74 Consequently, it is important to assess the ability of the inner AF to support cell injections, which cause a higher nutrient demand -and ultimately lead to low cellular viability.⁷⁵

Between Pfirrmann grades

Model predictions demonstrated that the degeneration grade strongly correlated with average nutrient concentration and critical solute levels (ie, glucose, oxygen, and lactate; Figure 5). As the disc degenerates, it undergoes morphological and biochemical changes.^12,13 IVD degeneration includes a significant loss in proteoglycans (PG), leading to disc dehydration, structural changes of the NP and AF, and calcification of the CEP.^12,13 Those factors were shown to inhibit solute diffusion in the disc causing glucose and pH to drop to their critical levels of 0.5 (mol/m³) and 6.7, respectively, making the IVD's environment harsh for cells survival.^17,24,40

The model showed that healthier discs (grades 2 and 3), with sufficient water content and intact ECM, could facilitate better diffusion of nutrients into the disc and improve clearance of lactate (Figure 5B and Table S1). Although the model suggested an overall trend of increasing lactate content and acidity in the disc as it degenerated, this was not observed as the disc degenerated from grade 4 to 5 (Figure 5B and Table 5). This may be attributed to the severe reduction in the height of grade 5 discs that possibly improved diffusion distance (Table 6) despite compromised solute diffusivity in the disc. Nonetheless, this observation does not indicate that a grade 5 disc is suitable for cell injections because a large portion of the disc's ECM has already been compromised.

Last, a large overlap exists between the quartiles of the same solutes in different degeneration grades demonstrating the limitations of assessing IVD nutrient status solely based on degeneration grades. For example, average glucose and oxygen concentration means in grades 3 and 5 were similar (Figure 5A1), but the distribution of said solutes was different such that the inner AF in grade 5 discs experienced deficient nutrient levels (Table S1). In addition, there is an overlap in average lactate values between grades 2 and 4, indicating that grade 2 IVDs might share similar lactate levels to grade 4 discs (Figure 5B1). This observation shows the importance of generating patient-specific nutrient distributions to accurately assess nutrient availability in the IVD.

Between patients

This study highlights the importance of making patient-specific model parameters such as disc size, metabolic rates, the diffusivity of the AF and CEP, and cell density. There was a significant correlation between increasing disc size to lower glucose and oxygen levels and higher lactate levels. For example, minimum glucose levels were 80% lower in a 50% larger IVD. This was expected since a larger IVD size means a longer diffusion pathway and vice versa.^75-77 A decrease in nutrient availability has been reported during the first decade of human life as the IVD increases in size.⁷⁸ Other factors such as surgery and degeneration can also alter the IVD's size leading to changes in nutritional status.

Model predictions show that degenerated discs with a relatively larger area exhibited deficient nutrients and high accumulation of lactate. For instance, IVD #18 in grade 4 has an area that is 27.5% larger than the group's average yet showcases minimum glucose levels below the critical concentrations at 0.39 (mM). Maximum lactate levels in this IVD are also high at 7.3 (mM). Despite this IVD not having any pathological conditions, its nutrient status cannot support normal cellular functions. In contrast, IVD #20 in grade 5 has a 45% smaller area than the group's mean, yet its minimum glucose levels are 2.5 (mM), and maximum lactate is 4.39 (mM), indicating nutrient conditions that can facilitate cellular functions. This improvement in IVD nutrient status may be attributed to the shorter diffusion pathway, which is also predicted in discs of varying sizes within Pfirrmann grades 2 and 3 (Figure 5B and Figure 7). In grade 2, IVD #2 and 8 have relatively small IVD areas compared to the group's mean and therefore have high nutrient concentrations. In contrast, patients #5 and 7 (grade 2) exhibited low average nutrient levels, yet their disc areas were not considered to be relatively large. In grade 3, only IVD #14 has a significantly large disc area, but surprisingly, it maintains a healthy nutrient status. Assessing the masks that were made for IVD #14 showed an artifact in the AF which did not enclose the NP completely which improved the nutrient status in this IVD. These observations point to the importance of tailoring the model's parameters to the specific IVD's geometry and size to more accurately assess nutrient status.

The SA shows glucose, lactate, and oxygen reaction rates significantly affect their corresponding concentrations in the IVD (Figure 8). Increasing the metabolic rates in the IVD leads to an increase in lactate levels and other by-products that are noxious to cells.^79,80 Aging and dehydration of the IVD, which can vary between patients, are two factors than can affect metabolic rates due to changes in cellularity and water content.⁸¹ Therefore, characterizing metabolic rates in patient IVDs with different degeneration grades may improve nutrient availability and solute transport modeling in the IVD. Different types of chemical exchange saturation transfer (CEST) MRI were developed to measure metabolism in vivo. CEST has also been optimized for human patients by significantly reducing the scan time to 5 minutes, rendering it viable to investigate metabolism in the IVD.^{82 31}P saturation transfer (ST) has been used to measure metabolic rates noninvasively in the human heart and brain.^83-85 This technique provides high spatial resolution and sensitivity in measuring energy metabolism in organs by assessing the conversion of phosphocreatine in creatine kinase reactions.^82-84,86 Magnetic resonance spectroscopy (MRS) has been applied to complement MRI techniques by providing information on the chemical changes that occur within the disc as it degenerates. For example, MRS has been demonstrated to quantify metabolites, specifically lactate and PG as markers for pain in the disc.⁷⁷ In a 2019 study, researchers were able to optimize MRS to accurately assess the integrity of the ECM by analyzing contents of PG, carbohydrates and acidity (alanine and lactate).⁷⁸ Incorporating these techniques into the modeling of solute transport have the potential in generating more accurate nutrient distributions in patient IVDs.

We also found a significant effect of AF and CEP diffusivities on solute concentrations in the IVD and the inner AF (Figure 8). We expected the AF's diffusivity to play a substantial role in determining critical nutrient levels because they lie in the inner AF. The effect of CEP diffusivity on the IVD's nutrient status has been addressed by several studies.^87-91 The CEP has a large number of channels and shares the largest surface area with the NP, which allows for a larger solute flux.^{29,67,87,90,92-98} In addition, the CEP is relatively thin (~2 mm) and porous, which facilitates better transport of molecules.⁹⁸ Due to the significant effects of AF and CEP diffusivities on solute concentrations in the inner AF and throughout the IVD, it is essential to include high-resolution ADC maps for both anatomical parts. Such maps have been generated to study variations in AF diffusivity in relation to degeneration grade and location within the AF (posterior, anterior, inner, and outer).^99-102

Interestingly, changes in NP diffusivity had a negligible impact on all measured metrics except for maximum lactate (Figure 8B). The NP facilitates better diffusion due to its biochemical structure and large water content.¹⁰⁰ Therefore, the chosen 50% change in NP diffusivity was not limiting to solute transport.

Cell density is shown to be a predictor for maximum lactate and average oxygen concentration (Figure 8B,C). Interestingly, cell density was not considered a predictor for glucose levels despite the strong correlation between cell density and maximum lactate (Figure 8B). We expected changes in cell density to substantially affect nutrient distributions in the IVD because it has been shown to alter metabolic rates.⁶⁴ It is suspected that the 150% dose of cells used in this model, which raises total NP cell density to 6000 cells/mm³, is too low to affect nutrient distribution compared to the highest dose used in some clinical trials which would increase total NP cell density to 11 000 cells/mm³ (Figure 10).^23,81 In addition, the study only assessed changes to the NP's cell density and overlooked any changes in the AF and CEP because cell injections often target the NP.^{23,81,103,104} It is also possible that because healthy discs were used for the SA that changes in cell densities did not have a large impact on IVD nutrient status. Further tests would need to be conducted to investigate the effects of changing cell density in the AF, and the CEP. Another untested but critical comparison is the effect of localized cell injection in the NP vs an increase in total NP cell density which could be useful in understanding the effects of cell delivery to various parts of the IVD.

FIGURE 10. Comparison of increasing the nucleus pulposus cell density on glucose distribution (mM) in patient # 2, grade 3. Cell density of 4000 cells/mm3 was used as a baseline value. 6000 cells/mm3 represents a 150% increase. While 11 000 cells/mm3 was the highest dose reported to be used by clinical trials23,82

Patient demographics showed IVD level to significantly vary between grade 2 and the remaining grades. However, disc level had no significant effect on minimum or maximum solute levels in the IVD (Kruskal Wallis, P <.05). The constant parameters chosen for this study, for example, water content, cell density, and so on, vary with IVD degeneration. Therefore, observed solute variation between groups could be attributed to degeneration grades. This makes sense given that healthier discs have shown to have better ADC values.²⁷ Lastly, nutrient distributions in patient outliers could not be predicted based on the degeneration level alone or IVD pathology, as was the case for patient #18. These IVDs had a specific disc shape and size, that led to specific nutrient distributions.

Last, the model predicted that modifications to disc size, metabolic rates, and the diffusivities of the AF and CEP greatly influenced the IVD's nutrient status (Figure 8). Therefore, future modeling techniques that aim to investigate the IVD's nutrient status should consider tailoring these parameters to the individual patient's IVD.

Model limitations

Some of the limitations to this model included possible artifacts in IVD masks. For example, IVD #14 in grade 3 did not exhibit abnormal levels of solutes in the disc due to a defect in the AF mask. This defect did not completely enclose the NP, and this significantly enhanced diffusion through the disc. This mask defect was a limitation of the model which inaccurately led to low levels of lactate and high levels of glucose and oxygen in discs with large areas.

Furthermore, the only model parameters made patient-specific were disc geometry and NP diffusion coefficients. There are many other parameters that could be expected to vary between patients that were not included in this model. A list of parameters and the basis of expected variation between patients was summarized in Table 7.

TABLE 7 Summary of current vs expected variations in solute transport parameters

Parameter	Current model specificity	Expected variation between patients
IVD Morphology	Patient specific	Expected to vary based on patient height, biochemical composition, gross defects, degeneration, mechanical loading, and level
AF solute diffusion coefficients	Pfirrmann specific	Expected to vary based on biochemical composition, location, mechanical loading, and gross defects24-26,67,68
CEP solute diffusion coefficients	Pfirrmann specific	Expected to vary based on biochemical composition, location, mechanical loading, and gross defects21-26
NP solute diffusion coefficients	Patient specific	Expected to vary based on biochemical composition, location and gross defects27-29,44
AF solute boundary conditions	Uniform	Expected to vary based on capillary density, and solute perfusion58
CEP solute boundary conditions	Uniform	Expected to vary based on thickness of the boney endplate, capillary density, number of marrow cavities and gross abnormalities69
Solute initial condition	Uniform	Expected to vary based on IVD physiological conditions
AF diffusivity	Pfirrmann specific	Expected to vary based on spinal morphology, biochemical composition, location, mechanical loading, and pathophysiology64,66,70
CEP diffusivity	Pfirrmann specific	Expected to vary based on spinal morphology, location, and endplate pathophysiology, including capillary density, lesions, and calcification38,64,69,70
NP diffusivity	Pfirrmann specific	Expected to vary based on spinal morphology, location, biochemical composition, and pathophysiology64,66,70
AF cell density	Uniform	Expected to vary based on glucose and oxygen availability, acidity, location, and patient age42,60
CEP cell density	Uniform	Expected to vary based on glucose and oxygen availability, acidity, location, and patient age42,60
NP cell density	Uniform	Expected to vary based on glucose and oxygen availability, acidity, location, and patient age42,60
Cell metabolism	Based on glucose/oxygen concentrations	Expected to vary more widely based on signaling factors, senescence, and patient age105
Water content	Pfirrmann specific	Expected to vary based on biochemical composition, location, mechanical loading, and pathophysiology64,66,70

Note: The specificity of each model parameter is given as well as the expected source of variation for each parameter. Patient specific factors are highlighted in gray.

Such limitations include applying constant homogeneous cell density between the different degeneration grades. Liebscher et al has shown that cell density varies even within the same anatomical part of the IVD; however, the same study suggested cell density in the AF, CEP, and NP does not correlate to disc degeneration caused by normal aging.¹⁰⁶ Still, this limitation significantly altered model results, as demonstrated by the SA. Cell density is expected to decrease with increasing degeneration grade, which would lead to reduced metabolic rates. Studies reported changes in the IVD's metabolic rates to alter glycolysis rates and change the levels of by-products including lactate that could lead to further degeneration of the IVD.^79-81

This model also assumed constant water content in the AF and CEP. However, studies demonstrated variations in water content based on location within the same tissue, which affected local solute diffusion.³⁰ Future models could eliminate this limitation by including a water content map for each part of the IVD. Antoniou et al, used different MR modalities (ie, T₁, T₂, T_1ρ, MTR, and ADC) to calculate water content in the NP and AF.¹⁰⁷ Applying this method in future iterations of the model should provide more accurate and specific model results. It was also important to note that a variation exists in published values for water content % and cell density (Table 8), and thus utilizing MR scans to personalize some of these parameters, for example, water content, could produce more accurate results. Also, the model did not consider other factors that can impact nutrient transport such as changes in load, charge density, and varying cellular metabolic rates.^28,111,112

TABLE 8 Literature values of water content and cell density showing various values reported and used for IVD modeling

			Water content (%)			Cell density × 103 (1/mm3)
Year	Author	Model	CEP	AF	NP	CEP	AF	NP
2016	Zhu[68]	3D model	—	70-85a	85	—	9	4
2010	A. Shirazi[24	3D model	60	66–73a	80	15	4–16	4
2007	Iatridis[108]		—	66	86	—	—	—
2005	Soukane[46]	3D model	—	—	85	—	9	4
1998	Bartels[47]	Human in-vivo	—	—	—	15	9	4
1988	Urban[109]	Human in-vitro	—	—	75-85a	—	—	—
1985	Kraemer[110]	Human in-vivo	—	73-80b	76-90b	—	—	—

^aThe high value for the inner AF and the low value for the outer AF.

^bA high value for healthy young individuals, while the low value for older individuals.

This model also simplified factors affecting solute perfusion to include calcification of the CEP only, which was modeled by decreasing CEP water content by 18%.⁴⁵ It has been shown that disc degeneration alters the number of capillaries contacting the endplates.¹¹³ Such changes in solute availability at the disc endplates could invalidate our assumption of constant boundary conditions across different Pfirrmann grades. Future work would seek to evaluate endplate capillary densities in-vivo.

Time of day when MR scans were acquired was not controlled. This limitation could considerably affect solute diffusion in the disc because different patient postures and activities during the day were shown to affect disc geometry and water content.¹¹⁴ For example, MR scans of IVDs taken in the morning showed higher water content and increased IVD height^114,115—two factors that were demonstrated by this model to significantly impact solute diffusion in the disc (Table 6). In addition, this model was created in 2D instead of 3D. This allowed for the collection of a large MRI data set, which would have been prohibitive with the longer scan time associated with more slices. Furthermore, it decreased model computation times. While the model showcased variations in solute distributions between patients and degeneration grades, the results were not validated experimentally. As with other models that have been developed, it was challenging to validate model results due to a knowledge gap regarding important parameters that were needed to model solute distribution in patients. The SA performed in this study had successfully identified an important set of parameters that needed to be obtained in order to validate the results experimentally. One of which was the metabolic rates for lactate in the IVD. As discussed earlier, changes to the metabolic rates of lactate was shown to alter solute distribution in the IVD, and it was correlated to maximum lactate and minimum glucose (Figure 8). Studies reported changes in metabolic rates, especially an increase in lactate production to create harsh IVD microenvironments that led to tissue degeneration.^79,80

Despite the model's limitations, the results showed general trends that were consistent across all patients which demonstrated consistent model performance. The specific trends of each patient were unique suggesting successful integration of patient data. Comparison of the specific trends for each solute across all patients showed distinct patterns that highlighted the importance of the individual morphology and physiological conditions of each disc, even among discs of the same Pfirrmann grade. Pending further development, and inclusion of more patient specific factors, this model may allow clinicians to account for the impact of increasing cell activity or density in a nutrient-starved environment. By iteratively increasing cell concentrations until a certain threshold has been reached the maximum capacity of a specific disc for new cells cell be determined. If levels were already low enough to indicate cell death, these patients could be excluded as candidates for cell injection-based therapies and advised to seek traditional treatments. The patients identified as poor candidates for cell therapies may be able to seek other novel therapies as they become available.

CONCLUSION

This work described a method for incorporating patient data into a nutrient transport model of the IVD. Results showed a decrease in glucose and oxygen concentrations accompanied with an increase in extreme lactate levels in the disc as it degenerates. Disc size, AF and CEP diffusivities, metabolic rates, and cell density had significant impacts on solute distributions in the IVD. The model also demonstrated distinct diffusion behavior between patients, even between discs of the same Pfirrmann grade. The importance of the distinct disc morphologies and physiological environments of each patient to the diffusion gradients in the disc was readily apparent. Patient-specific models could allow clinicians to further personalize treatments to the patient, for example, choosing the appropriate dose of cells to be injected into the IVD. This functionality indicated that patient-specific models could prove valuable in a clinical setting when predicting patient outcomes or treatment options.

ACKNOWLEDGMENT

The authors thank Dr. Aaron Fields for his valuable insights and helpful discussion. The authors also thank Lindsay Benage for her help in revising the article.

CONFLICT OF INTEREST

The authors have no conflict of interest to report.

AUTHOR CONTRIBUTIONS

Ward Shalash, Sonia R. Ahrens, and MBG were involved with research design, acquisition, and analysis of data. Liudmila A. Bardonova and Vadim A. Byvaltsev were involved with research design and data acquisition. Morgan B. Giers and Ward Shalash drafted and revised the document. All authors have read and approved of the final submitted manuscript.