An 83-year-old male who underwent right hip hemiarthroplasty due to a femoral neck fracture induced by a sideway fall was recruited. The femoral head was extracted during surgery and fixed with 4% paraformaldehyde for 48 h. The patient and family were notified and consented to contribute to the femoral head for research after hemiarthroplasty.

The femoral head was scanned using a high-resolution micro-CT scanner (SkyScan 1276, Bruker, United States) at 42 μm resolution with 8-bit gray-level values. The x-ray source voltage was set to 100 kV, and the current was 200 μA. A rotation step as low as 0.8° was selected to obtain the finest resolution. Based on the micro-CT data, 8-bit bitmaps (BMP) of every cross-section along the femoral neck axis were reconstructed using the bundled software Nrecon 1.7.1 (Bruker, United States). DataViewer 1.5.6 (Bruker, United States) software was employed to orient the cross-sectional images perpendicular to the main compressive trabecular bone in the femoral head. Finally, in CTAn1.16.8 (Bruker, United States) software, a global threshold of 50–255 was chosen to segment the trabecular bone on the BMP images, and a cube with dimensions of 3 × 3 × 3 mm was taken as the volume of interest (VOI) from the domain of the femoral head, as shown in Figures 1A, C.

The VOI cube was first used to investigate the anisotropy effect by applying NP in different directions. Further, to investigate the effect of individual trabecular structure remodeling on fluid dynamics, a cohort with a low global threshold limit from 80 to 50 at an interval of 10 was chosen to segment out trabeculae under the same VOI, and their bone volume fraction (BV/TV) was recorded (Figure 1E). VOI trabeculae demonstrate higher BV/TV with a lower global threshold limit setting; therefore, the effect of using a different lower threshold when segmenting trabecular bone can be assessed (Muller and Ruegsegger, 1996; Barak and Black, 2018). The VOI trabeculae were saved as “stl” files for subsequent model optimization.

Since the raw “stl” VOI trabeculae derived from CTAn software were full of noises (e.g., tiny spikes, unconnected elements, and unclosed pores) and were too rough for computational modeling, some polishing steps were undertaken using Geomagic studio 2013 software (GEOMAGIC Inc., United States). Optimization includes the “Mesh Doctor” (to detect and fix kind of flaws automatically) and “Remesh” (to acquire uniformly distributed two-dimensional elements). They were carefully operated to maintain the trabecular structure (Figure 1D). Simultaneously, a group of element sizes was set for the convergence test before the formal experiments. After optimization management, the VOI trabeculae were saved as “stl” files again, ready for the CFD model establishment.

In the SpaceClaim module of ANSYS 19.0 (Ansys Inc., United States), the polished VOI trabeculae were enclosed by a 3.4 × 3.4 × 3.4 mm fluid box to represent the marrow region. Then, the fluid-solid surface models were volumetric-meshed as tetrahedral elements using Hypermesh software (version 14.0; Altair Inc., United States) and assigned to the outlet and inlet planes. In each CFD simulation, the fluid box was assigned one outlet with NP and the other five planes were set to the inlet with zero pressure (0 P), simulating the bone marrow environment under a standard atmospheric pressure (Figures 2A, B). Although there are some differences between various NPWT systems concerning the filler material used (foam vs. non-adherent antimicrobial gauze) and the connecting suction catheter (integrated with a pressure sensor vs. flat drain), the intensity of NP ranges from −50 to −200 mmHg in most clinical practice (1 mmHg = 133 Pa) (Normandin et al., 2021), depending on the wound condition (Meloni et al., 2015; Agarwal et al., 2019). In the present study, NPs of −80, −120, −160, and −200 mmHg were applied along three principal directions (along the X-, Y-, and Z-axis, respectively) for anisotropy investigation. To study the effect of trabecular structural remodeling on fluid dynamics, an NP of −80 mmHg was chosen as a representative (Costa et al., 2020; Shiels et al., 2021).

The trabecular bone was considered rigid and had no slip wall, and the simulated marrow was assumed to be a steady-state, laminar, and incompressible Newtonian fluid with a density of 1.06 g/cm³ and viscosity of 0.4 Pa∙s (Rabiatul et al., 2021). The Navier-Stokes equation was adopted to define the flow behavior of viscous fluids. However, the effect of gravity was neglected. An iterative method was used to solve the equation for a steady flow, and convergence was identified when the relative tolerance was less than 0.001. Our preliminary study suggested that an element size of 0.10 mm of the fluid zone and trabeculae is sufficient to converge a steady fluid flow solution (Figure 2C).

As shown in Figure 2D, we defined the WD to enable post-processing of the CFD simulation. Specifically, a series of cross-sectional planes were placed at an interval of 0.4 mm to represent the WD from NP sources to the distal region. Cross-sectional planes were placed perpendicular to the NP direction. The resultant contents of each plane, including the pressure and shear stress on the bone and the velocity of marrow, were expressed as the median (25th–75th percentile) and reported as mean and standard deviation (‾x ± s) for regression (Wu et al., 2021).

The pressure, shear stress on the bone, and the velocity of marrow were determined using CFD simulation, we then attempted to verify the effect of different scales of NP on BMSCs gene expression, proliferation and osteogenic differentiation. By combining CFD results and cytological experiment, we could speculate how NPWT works within the bone tissue.

As shown in Figure 3, mouse BMSCs were seeded separately in a 24-well plate at a density of 1×10⁴ cells per well and cultured for 24 h. After that, the plates were placed in cell incubator with air inlet channel sealed and outlet channel linked to negative pressure device. The NP was controlled and adjusted using a vacuum pump (VSD Medical Technology Co., Ltd., customized, China) with a pressure sensor. The applied pressures of NP treatment were −80 mmHg, −120 mmHg, −160 mmHg, and −200 mmHg with a continuous pattern of 2 h per day, total 2 days for each plate.

After NP treatment, BMSCs were collected and harvested, and total RNA was prepared using the GenCatch TM Total RNA Extraction Kit (Epoch Life Sciences, 1,660,050) for RNA-sequence analysis following the manufacturer’s instructions (Illumina, Hiseq 2500, United States). Sequencing reads were mapped to the rat genome rn6 using HISAT2, and tag counts were summarized at the gene level using SAMtools, which allowed only one read per position per length. Differentially expressed genes (DEGs) were analyzed using edgeR. Gene set enrichment analysis, including Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG), was performed using DEGs with clusterProfiler. GO analysis was visualized using topGO, and Fisher’s exact test was used for statistical analysis. KEGG analysis was performed using clusterProfiler, and the hypergeometric test was used for statistical analysis.

After NP treatment, cells were continue to cultivate. At day 1, 3 and 7, the different medium was removed and replaced by 1 mL medium containing 100 μL cell counting kit-8 (CCK8) solution for 2 h incubation at 37 °C. After that, the optical density (OD) values of the cells in different medium were measured at 450 nm using an enzyme-linked immunosorbent assay plate reader (Multiskan MK3, Thermo Electron Corporation, United States).

BMSCs cells were seeded in 24-well culture plates as described in section 2.5. After 7 days of culture in different NP (2 × 10⁴ cells/well), the qualitative and quantitative ALP activity of BMSCs were characterized using BCIP/NBT alkaline phosphatase color development kit (Beyotime, China) and alkaline phosphatase kit (Nanjing Jiancheng, China), respectively.

For the qualitative ALP activity, the cells were fixed in 4% paraformaldehyde for 15 min and gently rinsed with PBS. The ALP activity staining was performed by a BCIP/NBT Alkaline Phosphatase color development kit according to the manufacturer’s instructions. Images were recorded with an Olympus MVX10 MacroView (Japan). For quantitative analyses of ALP activity, cells were incubated in RIPA lysis Buffer for 60 min. Before the analysis of ALP activity, total protein concentration (mg/mL) and ALP activity (units/mL) were measured following the manufacturer’s instructions with an enhanced BCA protein assay kit (Beyotime, China) and ALP reagent kit, respectively. ALP levels were normalized to the total protein content.

12 days after culturing under NP, the extracellular matrix mineralization of BMSCs were evaluated by the alizarin red S staining and calcium contents. The BMSCs were fixed in 4% paraformaldehyde for 15 min and gently rinsed with PBS. Finally, the sample were then incubated with 40 mM alizarin red S (ARS) staining solution (Sigma) for 30 min at room temperature. Images were recorded with an Olympus MVX10 MacroView (Japan). For quantitative analysis, calcium contents were measured using a calcium colorimetric assay kit (Nanjing Jiancheng, China) according to the manufacturer’s instructions. All samples were examined at less three times.

The pressure, shear stress upon the trabeculae, and flow velocity of the simulated marrow were analyzed. Specifically, these dependent variables were estimated as a function of the distance to NP sources (i.e., WD). Before regression analysis, we conducted a curve-fitting exam using SPSS software (version 20.0; IBM Inc., United States) and found a better-matched exponential function between WD and the dependent variables compared to other potential functions; therefore, a non-linear, least squares regression was used to determine the coefficients of an exponential function in the form: Y = A · e B x , where A and B are coefficients, e is the natural constant (Euler number), and Y is one of the three observed dependent variables estimated as a function of WD (x). The resulting coefficients with a p-value and correlation coefficient R ² are included.

Cytological experiment results are presented as the mean ± standard deviation (SD). Statistical comparisons were performed using a one-way analysis of variance followed by a Bonferroni test for multiple comparisons. Values of p < 0.05 were considered statistically significant.

As Figure 4 shows, the closer the NP source, the higher the pressure and shear stress on the trabeculae and the fluid velocity. Pressure and shear stress share a similar concentration area; however, the maximum shear stress tends to be located at tiny trabeculae. The fluid flow had significant anisotropy under NP in different directions; when NP was along the mechanical axis (Z-axis) of the trabeculae, the fluid flow was steady and laminar. A turbulent streamline can be found near the NP source when the NP is along the off-mechanical axis (non-Z-axis) of the trabeculae.

Pressure, shear stress, and velocity shifted down significantly from the NP source to 0.4 mm and 0.8 mm WD planes, especially when the NP was along the off-mechanical axis. The average median pressure, shear stress, and velocity at each WD for different NP are shown in Figure 5. All dependent variables at each WD plane increased with NP, and the increment of the first two WD planes was still more significant than that of the others. The changing tendency seems to be steadier when NP is along the mechanical axis than along the off-mechanical axis, indicating the influence of anisotropy of trabeculae on fluid dynamics. Although there was a 2-fold increase in NP (from 80 to 160 mmHg), there was no corresponding 2-fold increase in the dependent variables.

From Supplementary Tables S1–S4, non-linear regression demonstrated an exponential expression with a high confidence coefficient (R² > 0.96) between dependent variables and WD, regardless of the NP direction and scale. The pressure, shear stress on the trabeculae, and fluid velocity were well defined by the WD.

We used an NP of −80 mmHg as an example to investigate the effect of trabecular structure remodeling on fluid dynamics because the results derived from other scales of NP would be found to be similar to the −80 mmHg induced outcomes.

Figure 6 shows that the concentrations of the pressure and shear stress are located near the NP source. With trabecular remodeling, the concentration area increases accordingly, but the concentration depth does not change significantly. Streamlines become sparse and turbulent, particularly in regions close to the NP source.

As Figure 7 shows, pressure, shear stress, and velocity shifted down significantly from the NP source to 0.4 mm and 0.8 mm WD planes, especially when NP is along the off-mechanical axis, which is similar to the aforementioned results. Pressure, shear stress, and velocity decrease at most WD planes with structural remodeling, except at the 0.4 mm WD plane. This decreasing trend was more evident in the shear stress and velocity when the NP was along the off-mechanical axis. There was an approximately 2-fold increase in BV/TV (from 13.07% to 24.97%, corresponding to global threshold limit of 80 and 60); however, the decrease in the dependent variables was not 2-fold accordingly.

Nonlinear regression also demonstrated an exponential expression with a high confidence coefficient (R² > 0.97) between the dependent variables and WD. The confidence coefficient tended to decrease with structural remodeling. The pressure, shear stress, and velocity were better defined by the WD on the resorbed trabecular bone (Supplementary Tables S5–S8).

We performed RNA sequence analysis to verify the effect of different NP scales on BMSCs gene expression. Volcano analysis (Figure 8A) identified 2071 DEGs after −80 mmHg treatment, of which 1127 genes were upregulated, and 944 genes were downregulated. There were 1875 DEGs after −120 mmHg treatment, with 995 upregulated and 880 downregulated genes. For −160 mmHg treatment, volcano analysis identified 1886 DEGs, which included 1058 upregulated genes and 828 downregulated genes. In addition, 2043 DEGs were identified in the −200 mmHg group, of which 1113 genes were upregulated, and 930 genes were downregulated. From the Venn diagram (Figure 8B) of the −80, −120, −160, and −200 mmHg NP groups, more DEGs were observed in the −120 mmHg group than in the other groups.

GO analysis (Figure 9) revealed significant enrichment of genes related to biological processes, molecular functions, and cellular components. The enriched items related to the osteogenic effects of the NP treatment were further analyzed. The enriched genes were assigned to the following categories: “regulation of response to stimulus,” “cell proliferation,” “collagen-containing extracellular matrix,” “extracellular matrix,” “calcium ion binding,” “extracellular matrix structural constituent,” and “signaling receptor binding.” Consistently, the upregulated expressions of osteogenesis-promoting genes, such as Clca3a2 (chloride channel accessory 3A2), Angptl7 (angiopoietin-related protein 7), Adamts3 (a disintegrin and metalloproteinase with thrombospondin motifs 3), BMP binding endothelial regulator (Bmper), Bmpr1b, Dll1, Wnt4, Mapk3, and Mapk4, were identified in the −80 mmHg treatment group. In the −120 mmHg treatment group, the upregulated expressions of osteogenesis-promoting genes were Clca3a2, Angptl7, Adamts3, Adamts15, Bmper, Bmp6, Dll1, Wnt4, Wnt11, Wnt2b, Wnt9a, Mmp15, Stab2, Stap2, Mapk13, Notch3, Vegfd, Fgf2, and Mapk4. Osteogenesis-associated DEGs in the −160 mmHg treatment group included Clca3a2, Adamts3, Icam1, Mmp15, Bmper, Axin2, Adamts2, Dll1, Map315k4, Wnt4, Fgf2, and Notch3. In the −200 mmHg treatment group, the osteogenesis-associated DEGs included Clca3a2, Angptl7, Adamts3, Wnt2b, Icam1, Mmp15, Bmper, Axin2, Dll1, Heg1, Wnt4, Adamts13, Mapk13, and Wnt9a (Figure 10).

These results suggest that NP treatment systemically induced osteoblast differentiation in BMSCs, and the osteogenic effect was more evident in the −120 mmHg NP group. To explore the initial driver of the osteogenic effect of NP treatment, we performed KEGG enrichment analysis (Figure 11). We noticed that the “PI3K-Akt signaling pathway” as PI3K and Akt are the critical sensors of intracellular signal transduction pathways that promote metabolism, proliferation, cell survival, growth, and angiogenesis in response to extracellular signals. In PI3K-Akt signaling pathway, as showed in Figure 12, the DEGs enriched in the pathway activating molecules such as GF, ECM, ITGA, BCAP, PI3K, NUR77 and Bcl2, and thus activated PI3K-Akt signaling pathway.

Noticeably, in −120mmHg vs. Control group, the molecules of AMPK and CytokineR were significantly activated, which may cause to be different from the others groups. In addition, as Figure 13 shows, PLC, PI3K, p38, Bcl-2, VEGF, ICAM-1 molecules were upregulated in AGE-RAGE signaling pathway in diabetic complications; and as Figure 13B demonstrates, the Activin, PI3K, p38, Axin, Fgf2, Lefty2, IGF, Isl1 molecules were upregulated in signaling pathways regulating pluripotency of stem cells. We suppose that osteoblast differentiation in BMSCs treated with −120 mmHg may also be regulated through “AGE-RAGE signaling pathway in diabetic complications” and “signaling pathways regulating pluripotency of stem cells”, which were not presented in the others NP scale.

The cell viability of BMSCs increased significantly after NP treatment, especially within 3 days (Figures 14A, B). ALP is an important marker in the early stage of osteoblastic differentiation. ALP staining showed that compared with NP of −80, −160 and −200 mmHg groups, the −120 mmHg had more ALP-positive cells, which suggested that osteogenic differentiation was promoted (Figure 14C). As shown in Figure 14E, ALP activity was maximal in −120 mmHg (6.315 ± 0.485 IU/g) compared with the −80 mmHg (5.037 ± 1.353 IU/g), −160 mmHg (2.096 ± 0.879 IU/g) and −200 mmHg (1.975 ± 0.427 IU/g). Calcium nodule formation is typically used as a marker to evaluate the later stage of cell osteogenesis. Consistently, more calcium nodules were observed in −120 mmHg group after 12 d of NP treatment (Figure 14D). The quantitative analysis of calcium nodules was shown in Figure 14F.