extfieldsofstudy:
Materials Science
Medicine

s2fieldsofstudy:
Biology
Materials Science
Engineering

text:
A step toward engineering thick tissues: Distributing microfibers within 3D printed frames

Abstract Microfiber mats for tissue engineering scaffolds support cell growth, but are limited by poor cell infiltration and nutrient transport. Three‐dimensional printing, specifically fused deposition modeling (FDM), can rapidly produce customized constructs, but macroscopic porosity resulting from low resolution reduces cell seeding efficiency and prevents the formation of continuous cell networks. Here we describe the fabrication of hierarchical scaffolds that integrate a fibrous microenvironment with the open macropore structure of FDM. Biodegradable tyrosine‐derived polycarbonate microfibers were airbrushed iteratively between layers of 3D printed support structure following optimization. Confocal imaging showed layers of airbrushed fiber mats supported human dermal fibroblast growth and extracellular matrix development throughout the scaffold. When implanted subcutaneously, hierarchical scaffolds facilitated greater cell infiltration and tissue formation than airbrushed fiber mats. Fibronectin matrix assembled in vitro throughout the hierarchical scaffold survived decellularization and provided a hybrid substrate for recellularization with mesenchymal stromal cells. These results demonstrate that by combining FDM and airbrushing techniques we can engineer customizable hierarchical scaffolds for thick tissues that support increased cell growth and infiltration.

Although airbrushing is not as widely used as electrospinning, it offers significant benefits over the latter, especially in the context of translatable tissue engineering technologies and integration with existing manufacturing methods like 3D printing (Daristotle, Behrens, Sandler, & Kofinas, 2016). Airbrushed fiber mats are fabricated using a commercial airbrush that costs a fraction of the price of an electrospinning setup. This process mimics industrial-style fiber mat fabrication and is much faster, taking only seconds to minutes to build a mat of significant thickness with deposition rates of 0.1 ml/min or more Medeiros, Glenn, Klamczynski, Orts, & Mattoso, 2009;Oliveira et al., 2011). The significantly lower fabrication time reduces associated costs and lead time on device production. Fibers produced by this method can be deposited on any substrate, including in vivo, providing much greater flexibility in application over electrospinning (Behrens et al., 2014).
However, a challenge faced by both electrospun and airbrushed fiber mats is the ability to promote cell infiltration and nutrient transport in mats of a significant thickness (Rnjak-Kovacina & Weiss, 2011). Impaction of the fibers upon the deposition surface limits the overall porosity by compressing the structure into a nearly twodimensional sheet. This results in limited three-dimensional volume and smaller pore sizes while also being weak mechanically and difficult to handle. These are areas where the incorporation of 3D printing techniques is beneficial.
The flexibility and speed of airbrushing allow for nearly seamless integration with 3D printing technologies (Behrens et al., 2014;Chen, Townsend, Sell, & Martin, 2017). 3D printing has emerged as a means to rapidly create centimeter-scaled tissue scaffolds with customizable shape and structure (Chia, 2015;Hutmacher Dietmar et al., 2001).
Fused deposition modeling (FDM) is attractive for its relative ease of use, structure customization, rapid prototyping ability, low cost, and wide polymer compatibility, including biocompatible materials like polylactic acid (PLA), polycaprolactone (PCL), and polylactic-coglycolic acid (PLGA). However, the ratio of strut size to pore size in these scaffolds is inefficient for cell culturing which limits their promise in tissue regeneration applications. The minimum filament diameter producible from FDM machines is on the order of 200 μm while spatial precision in the x-y plane is near 10 μm. The layer thickness, or z resolution, is limited by the nozzle thickness and is on the order of 60 μm (Kuznetsov, Solonin, Urzhumtsev, Schilling, & Tavitov, 2018).
The optimal pore size range for tissue engineering is considered 50-300 μm (Karageorgiou & Kaplan, 2005;Oh, Park, Kim, & Lee, 2007). These pore sizes are theoretically achievable with FDM, however, the porosity in such a scaffold would be below 60% as constrained by the diameter of the printed struts. Large, open macroporous structure of 3D printed scaffolds provides minimal surface area for cell attachment and growth. Airbrushing can be used to provide additional surface area within the 3D printed scaffold macropores for cell growth and development while minimally affecting the overall porosity of the devices. The framework of the 3D printed struts separates the airbrushed layers while maintaining the porosity of the fiber mat. Furthermore, cells rely on signaling cues and attachments points from the underlying matrix which are absent in scaffolds printed by FDM but could be provided in airbrushed fiber mats. Overall, the benefits of and challenges faced by FDM make it a perfect pairing for airbrushed fiber mat integration (Chen et al., 2019;De Mori, Peña Fernández, Blunn, Tozzi, & Roldo, 2018;Mellor et al., 2017;Moroni et al., 2008;Rampichová et al., 2018).
Based on these considerations, we have developed a method of creating hierarchical scaffolds combining airbrushed microfiber mats periodically interspersed within a FDM 3D printed support as a way to open up the fiber mat structure to allow for improved cell penetration while maintaining the rapid manufacturing and customization and improving the cell interactivity of 3D printed scaffolds. The airbrushed mat is formed out of a biodegradable tyrosine polycarbonate that has been used for bone regeneration in the past and degrades on the time frame of months allowing for the natural ECM to take over (Jinku et al., 2015;Kim et al., 2015;Zhang et al., 2016). Culturing human dermal fibroblasts upon the scaffold prior to decellularization allowed for the functionalization of the hierarchical scaffold with cell-derived ECM to create a hybrid natural and synthetic structure. These hierarchical scaffolds show cell penetration, extensive fibronectin deposition from cells cultured on the scaffolds, durability during processing, and improved tissue integration in vivo. Using this approach, we have harnessed the benefits of both airbrushing and 3D printing methods to create a versatile scaffold for engineering thick tissues.

| Airbrushed fiber mats
Solutions for airbrushing were prepared by dissolving E1001 (1k), a tyrosine-derived polycarbonate, in 3 ml of tetrahydrofuran (THF) to get to a range of final concentrations from 5% to 10% w/v. In the notation Exxyy(nk), the xx and yy are the percent mole fractions of desaminotyrosol-tyrosine (DT) and poly(ethylene glycol) (PEG) with the molecular weight, n, in kDa; the remainder is desaminotyrosoltyrosine ethyl ester (DTE; Magno et al., 2010;Yu & Kohn, 1999). Solutions were mixed on a rotating shaker plate until all polymer dissolved.
The polymer airbrushing solution was loaded into the holding chamber of a commercial painting airbrush (Master Airbrush, G222-Set, 0.5 mm nozzle, dual-action, gravity fed). The airbrush was connected to 25 psi dry nitrogen source and the tip was stored in THF when not in use to prevent clogging. With the dual-action airbrush, depressing the trigger began airflow over the airbrush tip, while pulling back on the trigger initiated solution flow and fiber deposition. During optimization, airbrushing distance was varied from 10 to 20 cm, time from 3 to 9 s, and nitrogen pressure from 15 to 25 psi. The tensile modulus of the airbrushed and randomly oriented electrospun fiber mats was measured using a dynamic mechanical analyzer (Sintech 5/D, MTS Systems Corporation) set up for tensile testing. The molecular weight of E1001 (1k) before and after airbrushing was measured by gel permeation chromatography (GPC, Waters 515 HPLC Pump, 717 Autosampler, 2414 Detector, Agilent PLgel 5 μm 10 Å, 2 columns) relative to a polystyrene standard (n = 1).

| Scaffold fabrication
Orthogonal log-stack scaffolds for 3D printing were designed using CAD software (Sketchup, Trimble Inc.) with 1.5 mm spacing between each strut, with strut dimensions of 0.5 × 0.25 × 27 mm. Each scaffold was 10 layers tall with a 90 rotation between each layer to make a box pattern. Open source software (ReplicatorG) was used to convert the STL file to a G-code for printing. A pause command was added after every other layer to allow for airbrushing of the scaffold in situ. Scaffolds were printed using a commercial fused deposition

| Scanning electron microscopy
Sections of the scaffold were attached to aluminum sample pegs using double-sided conductive adhesive tape. Sample pegs were placed on a charge reducing sample holder, which was loaded into the SEM (Phenom ProX, Phenom World, 10 kV). Images were acquired at three different locations at low (255×) and high (1,000×) magnifications.

| Micro-computed tomography
Scaffold with and without airbrushed fiber mats were imaged on a Skyskan1172 at the Rutgers University Molecular Imaging Center (n = 1). Images were acquired at 34 kV, 210 μA, 0.3 step size, and voxel dimension of 1.92 μm. Images were reconstructed with CTAn (Bruker) and 3D rendered with FIJI.

| Fiber mat characterization
Coverage of the airbrushed fiber mats was determined from SEM images using ImageJ thresholding to differentiate fiber from pore. The covered area was calculated from the thresholded image. Fiber diameters were measured manually using ImageJ from SEM images. The number of fiber defects was counted using the ImageJ cell counter plugin. Micro-CT images were processed in FIJI to determine the porosity of the hierarchical scaffolds. Porosity was determined from thresholded image stacks representing one unit cell of hierarchical scaffold or pore.

| Scaffold mechanical analysis
Mechanical properties of hierarchical and empty 3D printed scaffolds were evaluated using a dynamic mechanical analyzer (Sintech 5/D, MTS Systems Corporation) by compression testing. Each 10-layer scaffold was cut to 5 mm × 5 mm in size prior to testing (n = 3). Scaffolds were compressed to 50% strain or to the 1000 N load cell limit, whichever was reached first, at a rate of 5 mm/min. Young's modulus was calculated from the linear portion of the stress-strain curve.

| AFM analysis
Mechanical properties of airbrushed fibers were assessed using AFM (Nanosurf AG, Switzerland) and ACLA-50 contact-mode cantilever with a stiffness of 48 N/m (AppNano). The scaffold used for these tests (n = 1) was 3D printed struts with sparse, airbrushed E1001 (1k) fibers across 2 mm gaps, cut into a small piece, fitted into a small petri dish using silly putty, and then mounted on a microscope under the AFM scan head. The force spectroscopy was performed using a distance of 50 μm so the full behavior of the cantilever acting on the fiber can be captured. The speed of the cantilever was 1 μm/s and the setpoint was 100-400 nN. The maximum calibrated deflection was 730 nm. Several data points were taken on multiple fibers across the 2 mm gap.

| Cell culture and decellularization
Human Dermal Fibroblasts (HDF; ATCC) under passage 10 were expanded in high glucose Dulbecco's Modified Eagle's Media (DMEM; Gibco) supplemented with 10% fetal bovine serum (Atalanta Biologicals) and 0.1% gentamycin (Sigma). Sterilized hierarchical scaffolds were placed in individual wells of a nontissue culture polystyrene 24-well plate and seeded at a density of 1 × 10 4 cells/mm 2 in 0.5 ml medium. Scaffolds were placed on a rocker for 1 hr at 37 C, after which 1.5 ml of media was added and the plate was incubated without shaking for 7 days. For decellularization, a published protocol (Mao and Schwarzbauer, 2005) was used except that the second incubation in lysis buffer was for 90 min, 37 C. Matrices were stored in PBS at 4 C for at least 48 hr before reseeding with cells.
Human mesenchymal stromal cells (MSCs, Texas A&M) under passage 5 were expanded in Minimal Essential Medium Eagle-Alpha Modification (α-MEM; Gibco), supplemented with 10% fetal bovine serum (Atlanta Biologicals) and 0.1% gentamycin (Sigma). Human mesenchymal stromal cells (MSCs, Texas A&M) were plated in basal medium on scaffolds with decellularized ECM at 2 × 10 4 cells/well in a 24-well nontissue culture-treated plate. Cells were rocked gently on an orbital shaker for 1 hr before static incubation for a total of 24 hr.

| Immunostaining and imaging
After incubation, scaffolds were fixed in 4% paraformaldehyde solution (Affymatrix/USB) with gentle rocking for 15 min at room temperature and washed three times with PBS. Scaffolds were blocked with 5% goat serum (Sigma) for 1 hr at room temperature followed by TrueBlack autofluorescence blocking (Biotium) for 60 s according to the supplied protocol. After three PBS washes, scaffolds were incubated with a polyclonal fibronectin antibody (Invitrogen, PA5-29578; 1:100 dilution) in 2% ovalbumin (Sigma) for 1 hr at room temperature with slight shaking on an orbital shaker. The scaffolds were again washed twice with PBS. A secondary staining solution containing Hoechst (AnaSpec; 1:500) for nuclei staining, actin stain with phalloidin 594 (Invitrogen) at 1:20 dilution or rhodamine-phalloidin (Thermo Fisher) at 1:50 dilution where indicated, and Alexa Fluor goat anti-rabbit (488 or 633 where indicated; Invitrogen; 1:500) in 2% ovalbumin was added and incubated in the dark at room temperature with slight shaking for 1 hr. After staining, the scaffolds were washed twice with PBS and stored at 4 C until imaging.
Scaffolds were imaged on a multiphoton TCS SP2 confocal microscope. Maximum intensity projections and z-stacks through the depth of the scaffold were collected. Scaffolds were bisected to image the cross-sectional distribution of cells and fibronectin with maximum intensity projections. Sprague-Dawley rats (12 weeks old, 350 g, Charles River) were anesthetized by isoflurane (2% v/v) and treated with buprenorphine (general analgesic, 0.075 mg/kg s.c.), bupivacaine (local anesthetic, 2 mg/kg s.c.), and Baytril (antibiotic, 5 mg/kg s.c.). The dorsal region was shaved and sterilized with three iterative betadine and isopropanol washes. Up to four 10 mm cutaneous incisions were made, 10 mm from the midline, and pockets were bluntly dissected to 10 × 10 mm 2 . Sterile implants (n = 4) or controls (n = 2) were placed within the pockets, closed by a 9-mm wound closure clip. Animals were euthanized by exposure to CO 2 , 14 days after implantation.

| In vivo analysis
Scaffolds and the overlying skin were collected and fixed in 10% formalin for 48 hr at room temperature.

| Histology
Fixed samples were embedded in paraffin and sectioned (Rutgers University Office of Research Commercialization Pathology Services).

| Statistical analysis
Significance was verified using the Student's t-test and expressed as the mean ± standard deviation for a minimum of three measurements.
Statistical significance was shown by * for p < .05.

| Optimization of airbrushed mats
Airbrushed fiber mat morphology was optimized by tuning air pressure, polymer concentration, spray duration, and collector plate distance. Pressures below 20 psi did not produce adequate shearing forces at the tip of the airbrush nozzle to form fibers, resulting in an increased number of beading defects. Increasing the pressure above 30 psi resulted in fiber whipping behavior and excessive force and damage to the deposited fiber mat. Therefore, subsequent experiments used a feed pressure of 25 psi. The concentration of E1001 (1k) polymer solutions in THF tested ranged from 5 to 10 w/v%. Fiber diameter increased with solution concentration, ranging from 1.9 ± 0.5 μm at 5 w/v% to 3.4 ± 1.2 μm at 10 w/v% ( Figure 1c). Defects, defined as polymer beading along the length of the fibers, were observed to varying degrees at all concentrations. The number of fiber defects increased with solution concentration, ranging from 0.6 ± 0.1 defects per fiber to 1.1 ± 0.1 defects per fiber, for 5 w/v% and 10 w/v%, respectively (Figure 1d). Above 10 w/v%, the solution was too viscous pass through the airbrush nozzle without clogging. From these results, we determined the optimal concentration for E1001 (1k) fiber formation from airbrushing was 6% w/v to minimize fiber defects and reduce fiber diameter. The molecular weight of the E1001 (1k) was at 206 kDa before and after airbrushing ( Figure S1).
Once polymer concentration was set, the airbrushing distance (10-20 cm) and duration (3-9 s) were investigated to determine optimal fiber mat distribution. SEM images were thresholded in ImageJ to calculate the fiber coverage area and plotted against time and distance (Figure 1e). At a distance of 10 cm, the density of polymer fibers was highest at all time points investigated. Increasing the airbrushing duration from 3 to 9 s at 10 cm decreased fiber mat porosity from 40% to 13%. However, 10 cm of distance did not allow for sufficient solvent evaporation, resulting in a wet mat with film formation at 6 and 9 s. Increasing the distance to 20 cm was ideal for consistent fiber formation, with coverage of 77, 62, and 48% for 3, 6, and 9 s, respectively. Distances beyond 20 cm resulted in fewer fibers landing on the substrate and reduced control over the deposition location and fiber density. Therefore, the optimal airbrushing conditions for 33 kDa E1001 (1k) of 6 w/v%, 25 psi, 20 cm displacement and 3 s duration were used in scaffold fabrication for subsequent experiments.

| Hierarchical scaffolds
In order to increase the porosity of the airbrushed fiber mat, a hierarchical scaffold approach was developed by iteratively airbrushing E1001 (1k) and 3D printing PLA or PCL struts. The 3D printed portions of the scaffolds were simple orthogonal log stack structures with 1.5 mm spacing between parallel struts. E1001 (1k) has a T g of 90 C so printing with PCL at 110 C facilitated fusion between the struts and E1001 (1k) without damaging the fiber mats. In contrast, printing with PLA at 215 C caused the E1001 (1k) fiber mats to melt and fall away from the struts.  Figure 3c shows the ratios of force to deflection at various points in the sparse fiber mat prepared for the AFM measurement. The forces were <2 μN and the maximum deflection was 5 μm. These forces are comparable to the adhesion forces of cells on a variety of substrates indicating the fibers will provide sufficient strength for cell motility (Potthoff et al., 2012(Potthoff et al., , 2014. Measurements made close to the struts as well as away from the struts were not statistically different indicating that proximity to the support does not make fibers more rigid. These measurements suggest that the fibers produced by airbrushing are flexible enough to give tactile feedback to the cells that are seeded on the fiber mat.

| Cell proliferation in hybrid scaffolds
One-week HDF cell cultures on hybrid E1001 (1k) scaffolds were carried out to observe cell organization on infiltration into sterilized scaffolds. Cells were observed across the surface of the hierarchical scaffold, spanning the pores between the FDM struts (Figure 4a,b).
After 7 days, full-thickness maximum intensity confocal imaging of scaffold cross-sections showed the presence of cells with a deposited fibronectin matrix at each scaffold layer (Figure 4c). The presence of robust fibronectin matrix indicated that the cells remained viable throughout the thickness of the scaffold.

| In vivo cellular infiltration into hierarchical scaffolds
The formation of integrated tissues within large defects in vivo necessitates the infiltration of cells and vasculature into scaffolding. F I G U R E 5 In vivo H&E staining indicating differential cellular infiltration into scaffolds. (a) Schematic of scaffold cross-section plane. H&E stained sections of hierarchical (b) and airbrushed only (c) from 14-day rat subcutaneous explants. Scale bar is 500 μm F I G U R E 6 HDF infiltration and assembly of a robust ECM before (a) and after (b) decellularization on Day 7. Scaffolds were recellularized with MSC for 24 hr (c). Scaffolds were fixed and stained for fibronectin (green), nuclei (blue) and actin (red, c only). Scale bars are 200 μm (a,c) and 100 μm (b) of the unsupported airbrushed mat was a dense polymer bundle with no cellular infiltration or vascularization to support tissue development.

| Scaffold decellularization and recellularization
The functionalizability and the robustness of the hierarchical scaffold were tested by running standardized decellularization and recellularization protocols on scaffold cultured with HDFs . After 7 days in vitro, scaffolds were decellularized resulting in the removal of the fibroblast cells from the matrix while retaining the base E1001 (1k) fiber mats with ECM on top (Figure 6a, b). Fibronectin fibril density and morphology was maintained throughout the decellularization process. Upon recellularization with mesenchymal stromal cells, the decellularized scaffolds promoted cellular attachment and spreading (Figure 6c). Therefore, the hierarchical scaffolds can withstand standard tissue engineering protocols including decellularization and can be recellularized with application-specific cell types (Mao, Hoffman, Wu, Goyal, & Kohn, 2017).

| DISCUSSION
Hierarchical scaffolds were successfully produced by iteratively printing polymer struts while airbrushing polymer fibers between printed layers to create a 1 mm thick fibrous structure with improved cell integration properties. This process provides additional surface area for tissue regeneration within the 3D printed scaffold with >97% porosity between struts and~70% porosity overall. Airbrushing proved effective at forming microfiber meshes that could be easily incorporated into the 3D printed frames. At only 3 s of airbrushing per layer, the additional fabrication time for the scaffolds was on the order of minutes which aligns well with the rapid fabrication associated with 3D printing. This provides a significant improvement over combined 3D printing and electrospinning setups which can take up to five times longer to fabricate a single fiber mat layer (Medeiros et al., 2009).
Developing airbrushing methods for new polymers requires balancing multiple factors including solution concentration, target distance, and air pressure. When the polymer concentration was low, below 5 wt%, the airbrushed solution failed to form fibers due to insufficient time during flight for adequate solvent evaporation. This resulted in the formation of a thin, continuous polymer film coating the entire device, the end of which was seen in the higher number of defects per fiber. At high polymer concentration, near 10 wt%, the increased solution viscosity combined with rapid solvent evaporation clogged the airbrush tip, leading to sputtering and larger diameter fibers. At the optimal concentration of 6 wt% E1001 (1k), polymer fibers were observed extending from the tip of the airbrush, flying through the air and forming a nonwoven mat on the target substrate.
The resulting fiber mat presented the fewest defects and kept the small fiber size essential for promoting cell infiltration. Reduction in fiber defects and control of fiber diameter are important for maintaining the homogeneity of the scaffold, thus ensuring that cell behavior does not vary due to the underlying substrate and that mechanical properties are consistent throughout. Other groups have investigated how different models of airbrushes and nozzle diameters can affect fiber quality and consistency (Tutak, Gelven, Markle, & Palmer, 2015). While airbrushing was carried out by hand in this study, robotic automation of the process to work in tandem with 3D printing is certainly attainable.
The hierarchical scaffold approach described here ameliorates the problems associated with poor cellular infiltration that has long been a challenge for microfiber and nanofiber implants. In our system, the 3D printed supports take on the role of opening up the microfiber mats while increasing the volume, controlling the shape and protecting the fibers against mechanical wear. Controlling the 3D printed design along with the density of airbrushing during scaffold fabrication as shown in Figure 1 allows for creation of scaffold with pore sizes within the optimum range. Struts produced by FDM provided shape, stability, and spacing to the airbrushed fiber mats while the microfibers filled otherwise empty pore spaces, as shown in Figure  The layered architecture can be exploited to engineer anisotropic tissues with gradients of composition, functionality, and morphology.
While the microfiber network and cell-derived ECM within a hierarchical scaffold is intricate, it is also robust. Specifically, upon decellularization, the fibronectin network was maintained along with the underlying airbrushed fiber network. Decellularized fiber mats have been explored as a way to increase the bioactivity of substrates . This presents an intriguing potential to create hybrid scaffolds with both artificial polymer structure and naturally derived ECM functionality. Ideally, the airbrushed fiber mat will be replaced over time by natural ECM as the E1001 (1k) degrades and is resorbed. Providing cell-derived ECM on the scaffold through recelldecell demonstrated method of functionalizing what would otherwise be a blank template with cell cues and attachment points prior to implantation. HDFs were used in our study because of their ability to rapidly-produce a robust matrix, but could be substituted with different cell types for other applications. For example, prior work by our group has shown that cell-derived ECM from chondrocytes and osteoblasts on electrospun fiber mats promote either chondrogenic or osteogenic differentiation of hMSCs (Mao et al., 2017). E1001 (1k) is a promising polymer for scaffold fabrication in bone and cartilage regeneration. It has been used previously by our group for bone regeneration as a salt-leached porous foam scaffold in nonloadbearing applications where it increased ALP, OCN, and mineralized calcium production in vitro and promoted osteoconductive bone growth in vivo (Jinku et al., 2015;Kim et al., 2011;Kim et al., 2015).
By harnessing the ability to precondition hierarchical scaffold with osteogenic or chondrogenic cell-derived ECM we could produce scaffold with specific bone or cartilage regenerative properties. Trabecular bone is highly porous, around 80%, with a compressive modulus of 0.22-10.44 MPa which we were able to approximate in our hierarchical scaffold, achieving 70% porosity and 26 MPa compressive modulus (Misch, Qu, & Bidez, 1999;Renders, Mulder, Van Ruijven, & Van Eijden, 2007). In cartilage tissue the natural stiffness is much lower, in the range of 0.1-2 MPa (Zhang, Hu, & Athanasiou, 2009). Chemical functionalization of E1001 (1k) is a potential future route of exploration for increasing the bioactivity of the hierarchical scaffold as well.
Unfortunately, E1001 (1k) does not have the melt properties required to be 3D printable which limits its options for scaffold manufacturing.
Optimizing an airbrushing protocol for E1001 (1k) and incorporating it into 3D printed frames of PLA or PCL provides a rapid alternative to porogen leaching, electrospinning, or other traditional fabrication methods. The flexibility of our hierarchical scaffold fabrication method makes it possible to mimic the properties of tissues, such as bone and cartilage while providing natural cell cues by altering polymer used, scaffold geometry, fiber density, and cell seeding type.

| CONCLUSION
Hierarchical scaffolds were fabricated for tissue engineering applications by combining the ability of fiber mats to support cells with the mechanical strength, open pore structure, and customizability of the 3D printed frame. Airbrushing process variables including polymer concentration, air pressure, spray distance, and duration were optimized for E1001 (1k) to obtain reduced defect microfiber mats. These fiber mats were easily and rapidly interspersed within the 3D scaffold during the printing process. Combining the two fabrication methods allows these scaffolds to be highly customizable by employing a wide variety of materials in nearly limitless geometries. Cells easily penetrate the distributed microfibers in vitro as a result of spacing provided by the 3D printed scaffold. In vivo testing showed improved cell penetration and tissue integration into hierarchical scaffold relative to controls. In vitro cultured HDFs deposited extensive fibronectin matrix throughout the airbrushed fibers. Scaffolds, airbrushed fibers, and the cell-derived ECM were robust enough to survive decellularization and recellularization, increasing the functionality and versatility of this tissue-engineering platform. By combining FDM and airbrushing techniques, we have rapidly engineered customizable scaffolds that support increased cell growth and infiltration.

ACKNOWLEDGMENTS
We thank Nanosurf (Dr. Saju Nettikadan) for loaning us the AFM, and Dr. Edward Nelson for helping us with the measurements. This study was supported by RESBIO, the "Resource for polymeric biomaterials" funded by the National Institute of Health, National Institute of Biomedical Imaging and Bioengineering Award Number P41EB001046.
The authors would like to thank Daniel Martin for assistance with confocal imaging.


extfieldsofstudy:
Computer Science
Medicine

s2fieldsofstudy: Biology Engineering

text:
Neuro-musculoskeletal flexible multibody simulation yields a framework for efficient bone failure risk assessment

Fragility fractures are a major socioeconomic problem. A non-invasive, computationally-efficient method for the identification of fracture risk scenarios under the representation of neuro-musculoskeletal dynamics does not exist. We introduce a computational workflow that integrates modally-reduced, quantitative CT-based finite-element models into neuro-musculoskeletal flexible multibody simulation (NfMBS) for early bone fracture risk assessment. Our workflow quantifies the bone strength via the osteogenic stresses and strains that arise due to the physiological-like loading of the bone under the representation of patient-specific neuro-musculoskeletal dynamics. This allows for non-invasive, computationally-efficient dynamic analysis over the enormous parameter space of fracture risk scenarios, while requiring only sparse clinical data. Experimental validation on a fresh human femur specimen together with femur strength computations that were consistent with literature findings provide confidence in the workflow: The simulation of an entire squat took only 38 s CPU-time. Owing to the loss (16% cortical, 33% trabecular) of bone mineral density (BMD), the strain measure that is associated with bone fracture increased by 31.4%; and yielded an elevated risk of a femoral hip fracture. Our novel workflow could offer clinicians with decision-making guidance by enabling the first combined in-silico analysis tool using NfMBS and BMD measurements for optimized bone fracture risk assessment.

It is of major importance to reliably and rapidly obtain information on the dynamic in-vivo bone stresses and strains as these are fundamental parameters in mechanostat theory that characterizes the interplay between dynamic mechanical loading and musculoskeletal health 1 : Dynamic bone stresses and strains control bone remodeling 2,3 are indicators of bone fractures 2,4,5 and hence allow for tailored and monitored clinical treatment 6 . Additionally, dynamic stress analyses are essential to evaluate the mechanical properties of endoprosthetic implants in-situ during dynamic activities to improve implant design, implant durability and in turn, patient outcome 7,8 .
In this regard, bone diseases demonstrate a massive public health problem, especially bone mineral density (BMD) loss that eventually leads to fragility fractures which are associated with serious consequences in terms of morbidity, mortality, and socioeconomic burden 9 . Of the 1.7 million fragility fractures in 2011 in the US 10 and 3.5 million fragility fractures in the E27 countries 11 , the osteoporosis-related femoral hip fracture in particular has become a major socioeconomic problem: The femoral hip fracture has a relative incidence of approximately 23%, making it the most common fragility fracture and has caused an annual treatment cost of $19 billion USD in 2004 in the US and €37 billion in 2010 in the E27. Because of the association of BMD loss with age, these numbers are expected to increase by roughly 25% until 2025 as the elderly population (>65 yrs.) continues to increase.
Targeting the early identification of adverse bone metabolism, direct measurements of stresses and strains in or on living individuals' bones are impossible, from an ethical point of view, due to their invasiveness. Therefore, they are restricted to a limited number of superficial bone sites within cadaver studies using strain gauges 1,[12][13][14][15] .
In this respect, finite-element analysis (FEA) has been used extensively to address musculoskeletal research questions by means of non-invasive computational analysis; e.g., to study the mechanical behavior of human bones 13,[15][16][17][18] , bone remodeling [19][20][21] , bone adaption 22,23 and/or to estimate the progression of osteoarthritis in normal and overweight subjects [24][25][26] . Thereby, it has become a state-of-the-art method to assign the Hounsfield units (HU) of computed tomography (CT) imaging data to the finite-element mesh to capture the patient-specific bone density distribution 14,17 . This so-called quantitative CT (QCT) FEA 4,27,28 has become a key tool in the assessment of the bone status, e.g., via predicting its mechanical bone strength [4][5][6]18,29,30 . For example, in assessing the risk of BMD loss-related bone fractures, it has been shown that a subject-specific FEA that utilizes QCT for material property identification, is beneficial over the current clinical standard that utilizes dual-energy X-ray absorptiometry-based (DXA) BMD measurements for diagnosis and medication effect assessment in osteoporotic patients 4,6,[28][29][30] .
In contrast to FEA, patient-specific neuro-musculoskeletal multibody simulation (NMBS) seeks to optimize clinical decision-making by reliably predicting numerous patient-specific quantities in terms of musculoskeletal health from only scarce, non-invasive medical imaging and gait lab data 31,32 . NMBS models are regularly employed to enable rapid and reliable dynamic analysis over an enormous parameter space, to gain insight into human locomotion [33][34][35][36][37] or to evaluate surgical treatment options 38 by providing methods that quantify forces in muscles, ligaments, and joints while avoiding invasive interventions. Another application of these NMBS is the extraction of boundary conditions for FEA 23 or the co-simulation of FEA and NMBS that can be used to investigate, for e.g., the effects of gait modifications 39,40 .
However, established clinical diagnosing standards including QCT-based FEA lack the capability to estimate stresses and strains that arise due to the dynamic multiaxial loading of the bone under varying loading conditions as usually present in daily living activities and specific to every individual 5 to ultimately allow for the functional analysis of the bone. Moreover, FEA is usually limited to static or quasi-static loading cases and involves an immense computational cost that renders it computationally impractical for dynamic analysis of the entire neuro-musculoskeletal (NMS) system 1,23 . On the other hand, NMBS assumes bones as rigid bodies which makes the estimation of the field parameters stresses and strains impossible. In an attempt to resolve these issues, model reduction techniques have been employed within biomechanical FEA 8,41 . In the context of clinically relevant load cases that involve only small linear-elastic deformations, modally reduced flexible superelements models 42 (SEMs) were derived from finite-element models (FEMs) of bone segments and were implemented into NMS multibody models to allow for strain computation: Hughes et al. 43 observed lower tibial strains in individuals undergoing repetitive physical activity, which might suggest a lower risk of stress fractures in those individuals, due to positive bone remodeling as a result of increased physical activity. Gervais et al. 7 reproduced the fracture site of an osteosynthetic plate at the lateral, distal femur in accordance with a case study and reported that only the dynamic computation of stresses yielded realistic results.
However, the implementation from SEMs into NMBS is still in its infancy, and several limitations need to be overcome to arrive at a NMBS with clinical significance. In fact, till date, an experimentally validated study investigating the translation of QCT-based FEA into the workflow of a subject-specific NMBS does not exist: It is unclear how subject-specific QCT-based FEA translates into a subject-specific NMBS in terms of the prediction performance of osteogenic field parameters and computational efficiency.
Despite these challenges, it would be desirable to enable the prediction of physiological, i.e. dynamic, bone stresses and strains during physiological tasks to monitor the bone health and provide a numerical indication for early tendencies towards abnormalities that may yield the risk of fracture. Moreover, the high prevalence of femoral hip fractures stresses the importance of monitoring dynamic stresses and strains especially in the living femur that due to the surrounding soft tissue is difficult if not impossible to analyze; thus, further complicates the validation and translation of such in-silico approaches into clinical practice. Therefore, the question of if and to what extent state-of-the-art NMBS can predict those biophysical stimuli for assessing the healthy and the osteoporotic femur under physiological-like conditions remains open.

Neuro-musculoskeletal Flexible Multibody Simulation for Patient-specific Bone Failure Risk Assessment
In addressing these limitations, the objective of the present study is to augment the clinicians' diagnosis and monitoring repertoire by enabling the reliable and rapid prediction of patient-specific stresses and strains, solely based on gait lab and medical imaging data. Hence, our work introduces a clinically motivated workflow, enabling neuro-musculoskeletal flexible multibody simulation (NfMBS) that features dynamic stress and strain computation ( Fig. 1) by incorporating QCT-FEAs into conventional, subject-specific NMBS.
Overview of the novel workflow. Our study investigates the in-silico prediction of BMD loss-related fracture risk in bones under the consideration of structural and functional patient specificity, i.e., of material and structural properties as well as of musculoskeletal architecture and motion dynamics of the individual NMS system.
We generated a specimen-specific FEM of a fresh human femur specimen from a male donor (age = 58 yrs) from QCT scans with an included BMD calibration phantom. Next, we derived a modally reduced flexible SEM from the very same FEM by deploying the Craig-Bampton method 42 . Both the domain-specific simulation models of the femur, i.e. the FEM and SEM, were validated by conducting static and dynamic biomechanical tests on the fresh human femur specimen after computer model reconstruction. Finally, two simulation studies were conducted: The first simulation study resembled the bone strength experiments done by Orwoll and coworkers 4 , in which QCT-based FEMs were deployed for the risk assessment of femoral hip fractures, due to BMD loss, within a cohort study. Then, native and BMD loss-affected SEMs were implemented into our previously presented NMBS Figure 1. Workflow for neuro-musculoskeletal flexible multibody simulation (NfMBS) featuring dynamic stress and strain computation. The workflow is based solely on data that is typically available in clinical settings and arrives at a neuro-musculoskeletal multibody model with flexible bone segments: In conventional NMBS, motion capture and force plate data are captured in the gait laboratory, e.g., motion marker trajectories and force plate measurements that are specific to a motion scenario. Subsequently, these data are fed to computational multibody models of the neuro musculoskeletal (NMS) system that are created from medical imaging data. Depending on the specific question at hand-kinematic or dynamic-in-silico predictions are performed, enabling non-invasive analysis of NMS quantities. In an attempt to extend these analyses toward fragility fracture risk assessment in bones, NfMBS furthermore enables the prediction of stresses and strains that arise from the dynamic loading during clinically relevant motion scenarios by reducing QCT-based finite element models to superelement models for their integration into NMBS. The illustrations marked with * were taken from Carbone et al. 55   www.nature.com/scientificreports www.nature.com/scientificreports/ 0.018 mm for the vertical femoral head displacement. The computation of the FEM resembling the static experimental trial took 19.56 min CPU-time.

Superelement model validation.
As it concerns the prediction performance of the SEM in comparison to the FEM, the difference in the vertical femoral head displacement was very small with an RMSE of 0.0843 mm, over the entire simulation time. Similarly, the von Mises stress distribution as predicted over the whole bone surface by the SEM is in great agreement with the validated FEM. The deviation in the von Mises stresses computed by the FEM and the SEM were negligibly small with an average RMSE of 0.0478 MPa (range 0.0223 MPa-0.0790 MPa) over the complete simulation time (Fig. 3).
Owing to the modal reduction in the degrees of freedom (DoF) from the FEM (DoF = 370'362) to the SEM (DoF = 32), the CPU-time for the dynamic experimental trial over 5 s in real time, decreased considerably from 328.68 min for the dynamic FEA to 0.15 min for the SEM. Although the modal reduction (t = 13.44 min) as well as the stress and strain recovery (t = 14.01 min) add up to the SEM generation, the SE analysis outperformed the dynamic FEA by the factor of 11.91 in terms of CPU time. Moreover, the SEM generation is required only once; after which the SEM can be incorporated into a variety of loading cases.

Evaluation of bone quality parameters.
To obtain a measure of the hip fracture risk, computational bone strength experiments were conducted following the protocol by Orwoll et al. 4 . Orwoll and coworkers calculated the bone strength (the load at which the fracture occurred) as the applied load from the resulting load-deformation curve at 4% deformation (=4000 µstrain) of the femoral head with respect to the greater trochanter (Fig. 4A,B). This technique was reported to provide excellent predictions of femoral bone strength in cadaver laboratory studies (n = 51, r 2 = 0.80) and an excellent correlation between femoral bone strength and osteoporosis-related hip fractures.
Moreover, we modeled osteoporotic bone as described by Little and coworkers 46 , which provided insight into the change of osteogenic parameters during dynamic squat motion for the fracture risk assessment after BMD-loss leading to osteoporotic conditions in the previously analyzed femur specimen. The osteoporotic bone model was parametrically derived from the above described FEM by reducing the Young's moduli of the cortical and trabecular bone elements to 84% and 67% (referred to as the osteopenic C 84 T 67 ) or to 68% and 34% (referred to as the osteoporotic C 68 T 34 ) in the cortical and trabecular bone compared to the originally QCT-based Young's moduli assignment, respectively. To verify the osteoporotic bone models, we conducted numerical bone strength analyses according to the protocol by Orwoll et al. 4 . The load-displacement curves were then generalized to stress-strain curves to allow for fracture risk assessment under more complex and varying boundary conditions (Fig. 4C).
The results of the numerical bone strength experiments exhibit a very good agreement in the load-displacement curves as predicted by the FEA and SEM up to the bone failure criterion of ε Orwoll = 4%. Moreover, the healthy, native bone model (F Orwoll-FEA = 6390 N, F Orwoll-SEM = 6612 N) as well as the parametrically modeled BMD loss-affected bone models C 84 T 67 (F Orwoll-FEA = 5480 N, F Orwoll-SEM = 5646 N) and C 68 T 34 (F Orwoll-FEA = 4095 N, F Orwoll-SEM = 4110 N) agree with the classification criterions, according to Orwoll et al. 4 , which reported the bone strength to be an average of µ = 5939 N ± 1919 N for the non-fracture cases and µ = 3782 N ± 1563 N for the fracture cases. As expected, the weakened bone exhibits lower bone strength, which results in higher strains at lower stresses.
Dynamic stresses and strains during two-leg squat. The SEMs of the femur were then implemented into a musculoskeletal multibody model (Fig. 5A). For the validation of the musculoskeletal multibody model, www.nature.com/scientificreports www.nature.com/scientificreports/ please refer to the Supplementary Information on the validation of the musculoskeletal multibody model (Fig. SI2). The obtained flexible musculoskeletal multibody model was used to simulate a realistic dynamic two-leg squat, for which the von Mises stresses and the strains as defined by Orwoll et al. 4 were calculated as an indicator of hip fracture risk (Fig. 5B).
For hip fracture risk assessment, we implemented the baseline (native C 100 T 100 ) and osteopenic (C 84 T 67 ) SEMs into our NfMBS and compared the in stress-strain curves defined in Fig. 4C as the measure of fracture risk. This enabled us to gain insight into the effect of BMD-loss on the fracture risk measures while including not only specimen-specific bone properties but also NMS dynamics.
The simulation of the entire squat motion cycle of 4.5 s in real time required 38 s of CPU-time. Since the simulated squat motion is an activity of daily living, the bone strains remained in a safe zone, i.e., the strains never exceed the 4% fracture threshold throughout the whole motion cycle, as expected. The osteopenic bone (red curve in Fig. 5), in which the stiffness is reduced to 84% in the cortical and 67% in the trabecular bone exhibits a higher strain and therefore, lower bone strength as compared to the native bone (green curve in Fig. 5). The strain measure that is associated with bone fracture, i.e. the strains in the femoral neck as defined by Orwoll et al. 4 , increased by 31.4% owing to osteoporosis-related BMD loss and hence, yielded an elevated risk of femoral hip fractures.

Discussion
Early, accurate, and efficient bone density assessment is mandatory to tailor and monitor the clinical treatment of, e.g., osteoporosis 5,6,47 . It has been shown that patient-specific FEA that utilizes QCT for material property identification is beneficial over DXA-based BMD measurement for diagnosis and medication effect assessment in osteoporotic patients 2,3,29,30 . However, established techniques, including QCT-based FEA, lack the capability to estimate stresses and strains that arise due to the dynamic multiaxial loading of the bone under varying loading conditions as usually present in daily living activities and specific to the individual NMS system.
To address these limitations, we introduced a novel, clinically-motivated workflow to enable neuro-musculoskeletal flexible multibody simulation (NfMBS) that features dynamic stress and strain computation ( Fig. 1) by incorporating modally reduced flexible SEMs into conventional, subject-specific NMBS to enhance clinicians' diagnostic repertoire on BMD loss-related bone fractures. We present a validated workflow The load-displacement curves up to failure (ε Orwoll = 4%) are depicted for the native bone (green marker) and for the weakened bones (orange and red markers). Both, the FEA predictions (green curve) and SEM predictions (blue curve) are shown for validation purposes. For reference, the relative incidence of bone fractures (cyan for non-fracture cases, magenta for fracture cases) is depicted in dependence on the absolute bone strength as measured by Orwoll and coworkers. (C) The load-displacement curves are generalized to the material characterizing stress-strain curves as depicted for the superior, inferior, anterior, and posterior site of the femoral neck.
that establishes a general method which can estimate dynamic femoral bone stresses and strains, non-invasively. This method can help in developing bone mechanoregulation theories and in identifying various risk scenarios of femoral bone failure and specifically, femoral hip fractures.
The major contributions of our work are the integration of a subject-specific, QCT-based FEA of the human femur into the workflow of patient-specific NMBS together with its validation. The inclusion of structural patient specificity, i.e. morphological and material properties of the bone, and functional patient specificity, i.e. musculoskeletal architecture and motion dynamics of the NMS system, is necessary to enable clinically relevant NfMBS to identify early tendencies towards femoral hip fractures as it is a major risk scenario. We could show that our QCT-based FEM translated very well to the flexible SEM for its implementation into NfMBS, while maintaining the capability to reproduce osteoporosis-related bone weakening under dynamic NMS boundary conditions. The comparison of FEM and SEM exhibited negligible deviations between the target parameters at the femoral neck in terms of von Mises stresses and bone strength. In this context, our experimentally validated in-silico analysis of the human femur under loading conditions normally seen in daily living activities is the first study of its kind and provides great insight into highly prevalent femoral hip fractures, which are practically impossible to analyze in-vivo. Although the bone quality was varied only parametrically from healthy (C 100 T 100 ) to osteopenic (C 84 T 67 ) to osteoporotic bones (C 68 T 34 ), the results seem representative of osteoporosis-related fractures, since the predicted femoral bone strengths were in very good agreement with the literature data 4,5,48 . Our results provide confidence that NfMBS could predict the individual bone strength (as measured in Orwoll et al. 4 ) as a major diagnostic parameter for the assessment of risks and medication effects under physiological-like conditions in an accurate and efficient manner. Notably, is has been shown that there can be more complex criteria derived from the predicted mechanical stresses and strains for quantifying the resistance of the femoral bone against failure [49][50][51] .
Our most important limitation was that our study has a rather methodological character with its emphasis on the validity of the proposed workflow. A limitation of the presented work is that only one healthy femur specimen was analyzed. From the clinical point of view, it would be, therefore, desirable to apply the workflow to a larger, statistically relevant number of osteoporotic bone specimens and eventually to patient cohorts.
However, an extensive body of studies has shown that QCT-based FEA is a powerful diagnosis tool that possesses superior accuracy and reliability in the assessment of bone strength and bone fracture localization over DXA-based BMD measurements 5,6,29,30,47 . As mentioned above, BMD alone is not sufficient to describe the in-vivo intricacies that involve bone morphology and material-structural bone properties 47 . Thus, the logic extension to QCT-based FEA is the integration of patient-specific NMS dynamics. Clearly, the risk of bone fractures could be more accurately diagnosed with the complementation of NfMBS-predicted bone strength data under physiological-like representation with skeletal dynamics and soft tissue interaction 7,8 . NfMBS provides much deeper insight into relevant risk factors, as NMS dynamics-such as soft tissue damage or motor www.nature.com/scientificreports www.nature.com/scientificreports/ control dysfunctions-can be included into the set of boundary conditions to ultimately better resemble the patient-specific in-vivo situs. Consequently, NfMBS could enhance key aspects of QCT-FEA to ultimately arrive at a predictive tool that enables a combined in-silico analysis protocol of BMD measurement and dynamic NMS computer simulation for an improved diagnosis, e.g., of osteoporosis.
In fact, the direct comparison of our experimentally validated QCT-based FEM against the flexible SEM, in terms of the prediction accuracy and efficiency of clinically relevant bone parameters, revealed that NfMBS can handle the trade-off between these aspects extremely well, i.e., the SEM preserved very high accuracy while outperforming the dynamic QCT-based FEA in terms of magnitude in CPU-time. It should be noted that, currently, only NfMBS possesses the capability to sample the enormous parameters space for clinically relevant solutions by running numerous non-invasive parameter analyses on the patient's NMS system under dynamic boundary conditions 32 . Contrarily, very complex load cases, like the forward dynamic NMS squat simulation, have not been achieved up till date, by the use of FEA.
Regarding the QCT-based FEA, our FEM resembled the structure-mechanical properties of the tested bone specimen by, first, reconstructing the bone specimen's morphology using high-resolution QCT and then, discretizing the 3D solid into quadratic tetrahedral finite elements that were assigned bone-specific BMD-dependent Young's moduli for inhomogeneous and HU-value dependent, but isotropic and linear-elastic, material properties. The structural behavior of the bone in the plastic strain range beyond the yield point, especially the bone's failure mechanism itself, has not been investigated within this study and is not sufficiently described using a linear-elastic material model. However, the bone strength assessment according to Orwoll et al. 4 is applicable: the stress-strain curve of the femoral bone is roughly linear until the ultimate load is reached 17,18,48,52,53 , which should never occur during daily living activities. Generally, the consideration of microstructural processes along a timely trajectory, e.g. the architectural change of the trabecular bone, cannot be easily predicted using NfMBS. After capturing the patients current NMS condition, the clinically-relevant predictive capacity of NfMBS lies more in the incorporation of structure-mechanical and NMS boundary conditions for rapid and extensive osteogenic parameter analysis than in the inclusion of overly complex FEMs, as simple macroscopic QCT-based FEAs have been proven efficient 17,18,29,52,53 . Moreover, it was shown that there is only a marginal difference in the mechanical response when using anisotropic material models 54 .
However, our results show that decreases in BMD, e.g. due to osteoporosis or lytic bone changes, lead to abnormalities in the stresses and strains, which indicates a higher risk of fractures especially around the femoral neck. This agrees with numerous studies on the fracture risk of osteoporotic bone and its clinical prevalence, which report femoral hip fractures as the most prevalent osteoporosis-related fracture 4,9,48 . In this context, NfMBS resolves the trade-off between the predictive capability of fracture risk/fracture site and computational efficiency.
As such, our approach might be particularly valuable in the regular assessment of patient-specific bone strength, i.e., as an early diagnosis tool of fracture risk before the possibility of a fracture arises. In this respect, Anitha et al. 47 has shown that even low-quality CT imaging does not affect the prediction capability of QCT-based FEA. Hence, more frequent follow-up CT scans in combination with NfMBS would allow more regular monitoring and earlier treatment of BMD loss-related bone conditions.
Further limitations apply to the incorporation of flexible SEMs of bone structures into NMBS and the NMBS itself as follows: idealizations and simplifications of the real biomechanical system are inherent to NMBS that uses the multibody formalism to approximate the system's dynamics by means of rigid bodies, geometric constraints and discrete force elements. An overview can be found in 31,32 . We took great care of the representation of the musculoskeletal geometry, i.e. bone morphology, kinematic topology, and muscle attachment sites, throughout the entire modeling process, since musculoskeletal architecture is known to have considerable impact on muscle-driven motion patterns and joint loading 55,56 , and most likely on bone loading as well. For validation purposes, the entire workflow was specifically based on experimental datasets 45,55 , and the musculoskeletal architecture, in particular, was refined in close collaboration with an orthopaedic surgeon. Additionally, we validated the boundary conditions due to musculoskeletal dynamics in terms of the quadriceps force 57,58 and the resultant hip reaction force 59 (Fig. SI2).
As it concerns the neuro-muscular control, the computed muscle control algorithm may allow simulation of muscle-induced motion; however, the solution of using static optimization techniques for muscle distribution problems, although common in musculoskeletal multibody dynamics, is still a matter of debate in terms of neurophysiological evidence. In order to obtain dynamic stress and strain distributions for patient-specific and load case dependent scenarios, it is not only necessary to include bone-specific material properties but also subject-specific musculoskeletal architecture, i.e., muscle forces that act on the skeletal system as they would in-vivo 8 . Accomplishing this is still challenging as direct muscle force measurement is not possible and although statistical methods exist, automated implementation of musculoskeletal architecture is an unresolved problem 40 . Passive forces due to muscular or capsular structures were not considered, since passive soft tissue accounts for less than 10% of intersegmental torques during comparable activities 60 .
The use of publicly available datasets (SimTK 45 , TLEM 2.0 55 , Bergmann et al. 59 ) and the complementation of these data by our own experiments allowed for efficient validation of our workflow. This, however, stresses the importance of research on methods that help to identify an individuals' musculoskeletal system on the way toward clinically-relevant NfMBS 40 . Despite limited technical or methodological challenges, a considerable effort must be directed toward the clinical integration of QCT-based FE modeling and, ultimately, of NfMBS.
By defining a symmetry condition between the pelvis and the anatomical sagittal plane, we simplified the motion of the contralateral side. Kinematic studies suggest that this simplification is valid as humans show equivalent kinematics in the contralateral side during sitting down and standing up activities 61 . The two-leg squat, in particular, is a rather symmetrical motion and the symmetry condition in the sagittal plane showed only minor reaction forces in the constraint directions. This suggests that there was no abnormal bending of the femur.
www.nature.com/scientificreports www.nature.com/scientificreports/ In summary, we introduced a novel workflow for subject-specific NfMBS featuring dynamic stress and strain computation. The workflow aims at the efficient integration of the capabilities of QCT-based FEA into established subject-specific NMBS. Our work provides the foundation for a clinically motivated workflow that incorporates QCT-based FEMs into subject-specific NMBS. The performance of the novel workflow is supported by experimental results on a fresh human femur specimen and is in good agreement with comparable studies. Our results provide confidence that NfMBS can provide a computational tool of predictive in-silico analysis that could help in developing bone mechanoregulation theories, in identifying various risk scenarios and ultimately in optimal treatment options for femoral bone failure; especially in terms of highly problematic femoral hip fractures due to BMD loss. Until clinical translation, however, numerous challenges that are mainly associated with musculoskeletal modeling, e.g. an automatized implementation of subject-specific musculoskeletal architecture, need to be overcome. The application of NfMBS to a larger number of specimens, including osteoporotic bone specimens, is currently in progress. NfMBS seems promising in providing a generalized and reliable framework for the assessment of fracture risk and treatment options in osteoporotic patients.

Materials and Methods
For our NfMBS that features dynamic stress and strain computation, we used the TCP/IP communication framework that comprised the server SIMPACK ® (v9.7, Dassault Systèmes Deutschland GmbH, Gilching, Germany), in which the NMS simulation model with a flexible femur SEM was implemented, and the client Matlab/Simulink ® (v8.1, 2013a, The MathWorks Inc., Natick, USA), in which the computed muscle control was implemented. All simulations were run on an off-the-shelf desktop PC (Intel ® Core ™ i7-4790 CPU @3.60 GHz 8.00 GB RAM).
Experimental testing of a fresh human femur. We conducted static and dynamic axial compression tests on a thawed fresh femur with the experimental protocol adapted from 15 (Fig. 6). The specimen was received through an anatomic gift program (Science Care Inc., Phoenix, Arizona, USA) and was authorized by the ethical review committee (Ethics Committee of the Bavarian Medical Association, BLAEK 2011-058, Munich, Germany). We confirm that all experiments and research were performed in accordance with relevant guidelines and regulations, and informed consent for study participation was obtained from the donor. The experimental protocols for measurements on the human femur specimen were approved by the cluster 'Numerical Simulation' within the Musculoskeletal Biomechanics Network (MSB-Net) of the Basic Research section of the German Society for Orthopaedics and Trauma (DGOU).
We measured the global deformation, i.e. the displacement of the force application point, and the local strains, i.e. the strain distribution in the cortical bone, in response to the compression within the bone's elastic deformation range by using linear strain gauges. The detailed experimental protocol can be found in Supplementary Experimental Protocol of the in-vitro analysis. Note that mechanical testing is not mandatory to the proposed workflow and served validation purposes only.
Quantitative CT scan of a fresh human femur specimen. As depicted in Fig. 6, the femur, with its entire soft tissue structures removed, underwent a quantitative CT (QCT) scan for image-based FEM and SEM generation. The CT scan was conducted at the Trauma Center Murnau (SOMATOM Definition AS + CT scanner, Siemens AG, Erlangen, Germany) and included a BMD calibration phantom (European Forearm Phantom-05-83, QRM GmbH, Möhrendorf, Germany) for image-based material property assignment and calibration purposes. Specifically, the phantom served as quantification of the distributed bone stiffness via relating material specific X-ray attenuation in Hounsfield units (HU) to BMD and ultimately to the distributed Youngs' moduli. www.nature.com/scientificreports www.nature.com/scientificreports/ The reconstruction of the CT scans to the 3D FEM followed the protocol described in 14 . The CT images were saved as DICOM files with a resolution of 0.29 mm × 0.29 mm and a slice thickness of 0.6 mm, resulting in 886 layers. The DICOM files were imported in the segmentation software package AMIRA ® (v5.4.1, FEI Visualization Sciences Group, Hillsboro, Oregon, USA) to reconstruct the 3D bone surface model using triangulated surfaces. The resulting surface model was smoothed and converted into a surface of analytical non-uniform rational B-splines using Geomagic Studio (v2013, 3D Systems, South Carolina, U.S.A.), and finally exported to the FE pre-processor via the IGES interface.
As we were considering clinically relevant load cases with rather small deformations, we chose a linear-elastic HU-dependent material model. More specifically, our in-house mapping algorithm assigns a Youngs' modulus E i to each node i of the FEM that is in dependence on the HU values obtained from the QCT scans of the very same femur at the respective node location. To do so, the calibration phantom with a known BMD served as a reference for the linear regression relating the HU value of the phantom to the known BMD's equivalent ρ ash , which is the so-called ash density: .
The bone's Young's modulus could then be assigned nodewise, according to its respective ash density, following another relation presented in Cong et al. 62 : where E i is the Young's modulus in GPa and ρ ash,i is the ash density in g/cm 3 for the respective FE-node i. A Poisson's ratio of 0 3 ν = . was assumed over the entire bone 46 . Note that values below 1 MPa were set to 1 MPa, to avoid negative or zero elasticities 63 .
The mesh quality was optimized within a convergence analysis: the element size was incrementally decreased from 7.4 mm to 2.0 mm until the change in the computed strains and deformations fell below 5%. The final mesh for our baseline model consisted of 123,454 nodes and of 83,845 elements with an element size of 3.8 mm.
The numerical boundary conditions and locations of the strain gauges were extracted from the 3D scan to ensure that they resembled the experimental setup. To go into further detail: The FE nodes of the distal part of the femoral cortex below the casting resin surface were entirely constrained to zero DoF, and a compression load similar to the experimental loading was applied to a reference point that kinematically distributed the resulting displacement to the FE nodes on the upper hemisphere of the femoral head.
Based on the described static FE model, the femur specimen was also analyzed within a dynamic FEA with an implicit solver over 10 s simulation time. A dynamic loading regimen was defined and a constant mass density of 1940 kg/m³ was assumed 64 . The dynamic loading regimen was implemented as sinusoidal function equivalent to the experimental loading regimen as described in Supplementary Information on the experimental protocol of the in-vitro analysis. The resulting force was applied at the femoral head to a reference point that was directly coupled to the FE nodes of the femoral head via kinematic coupling representing the experimental stamp. Superelement model. Since clinically relevant load cases only cause small linear-elastic deformations in the bone that otherwise undergo large rigid body motion, the Craig-Bampton modal reduction method 42 in the floating frame of reference formulation was employed for the SEM generation from the specimen-specific FEM.
The unreduced FEM, is converted into the SEM by projecting the large number of FE-related DoFs into reduced matrices of mass ∼ M , stiffness K ∼ and modal coordinates u  using Craig-Bampton modes, The mode shapes CB Φ account for static and dynamic deformations due to forces, constraints, and joints acting upon the SEM, but this is, however, only within a limited frequency range. As a result, the mode shapes Φ k CB, that correspond to higher frequency responses of the SEM can be truncated without a relevant loss of information. The reduced SEM, i.e., the matrices M , K ∼ ∼ , the modal coordinates  u, and the modes shapes Φ CB , were calculated within a Craig-Bampton mode shape analysis using Abaqus/CAE ® .
Finally, the flexible bone structure was integrated via the flexible body input interface into the multibody dynamics software environment SIMPACK ® for subsequent, flexible multibody dynamic analysis. The numerical boundary conditions were defined similarly to the FEM (Fig. 7). All interface nodes (i.e. the respective reference point in the FE model) were connected to the respective subset of coupling nodes (i.e. the respective nodes of the FE model for the definition of the respective boundary conditions) using structural-continuum distributed coupling, which allows the coupling nodes to move relative to each other and reduces the impact of boundary conditions on the calculation of stress and strain. www.nature.com/scientificreports www.nature.com/scientificreports/ As required for dynamic multibody analysis, a constant mass density of 1940 kg/m 3 equivalent to the dynamic FEA was assumed 64 , the number of mode shapes was set to 30, and a proportional damping ratio of 0.025 was selected for all modes. Note that the number of mode shapes and the proportional damping ratio were verified by means of the dynamic compression tests and bone strength experiments.
As it concerns the stress and strain computation, the SEM's stresses are recovered by deploying a linear-elastic material law. According to 65 , the (6, 1) vector of stresses σ in a point of the flexible body is related to the (6, 1) vector of strains ε by σ ε = E with the symmetric (6, 6) matrix of elastic coefficients E. With the strain vector ε defined by the partial derivatives of the deformation vector u according to ε = Du using the (6, 3)-partial derivative operator matrix D and Φ = u u CB  with the Craig-Bampton shape functions CB Φ used in (4), the stress vector σ is related to the modal coordinates  u by Musculoskeletal multibody model. The NfMBS model (Fig. 8) was based on two complementary experimental datasets. The first dataset consisted of the experimentally validated SEM of the human femur specimen, as described above. The second dataset was obtained from the SimTK website (https://simtk.org/home/kneeloads) as made available to the public within the 4 th Grand Challenge Competition to Predict In Vivo Knee Loads 45 . To put it briefly, the SimTK dataset is a comprehensive and consistent collection of subject-specific models and measurements from one male subject (age = 88 yrs, height = 168 cm, m = 66.7 kg) that enables the validation of NMBS estimates in the lower extremity. The database included the geometry of the lower-right extremity (pelvis, femur, patella, tibia, fibula, and pes) from CT scans (pre-and post-op) and marker trajectories from motion capture experiments. To implement the femur specimen SEM into the NMBS, bone geometries and marker trajectories of the SimTK dataset were globally scaled to the size of the femur specimen (111.53%) using Geomagic Studio.
The reconstructed bone segments were mutually linked by ideal joints to an open kinematic chain (Fig. 8A). Starting from the ground, the inertia-fixed pes was linked via a two DoF universal joint to the tibia and fibula, which were summarized as one single segment. Similarly, the tibia-fibula compound was connected via one DoF hinge joint to the flexible femur SEM, which was in turn connected to the pelvis via a three DoF spherical joint. Finally, the upper body 66 was linked via a one DoF revolute joint to the pelvis' sacrum endplate center. In addition, spring-damper force elements (referred to as symmetry condition) restricted the motion of the sagittal plane by exerting forces and torques to the pelvis with respect to the ground, thus, eliminating abnormal bending torques around the sagittal plane during a symmetric squat motion. In this manner, the lower-left extremity was represented by its reaction forces and torques, acting upon the pelvis.  Fig. 6. The finite-element model (A) is modally reduced using shape functions obtained from a Craig-Bampton 42 mode shape analysis (B) to a flexible superelement model (C). In the finite-element model, the clamping is modeled by a zero degree of freedom joint. Implemented into the multibody simulation framework, the superelement model interfaces the multibody system via interface nodes on which forces and torques are exerted upon. In the model of the experimental setup shown in (C), the interface nodes are the load force application point 1 and the fixed clamping 2. In the musculoskeletal model shown in Fig. 8, interface nodes are provided for each joint marker, each muscle attachment site, and each via-point, respectively. www.nature.com/scientificreports www.nature.com/scientificreports/ The preservation of the bone geometry and anatomical landmarks allowed for the definition of mass properties, for the identification of joint rotation axes, by fitting spheres or cylinders into the articulating surfaces of adjacent segments, and moreover, for the establishment of standardized reference systems 67,68 , which allowed for the description of joint dynamics and the identification of the attachment sites of relevant muscle structures of the lower extremity as documented in the Twente Lower Extremity Model 2.0 dataset 55 . The origins and insertions of the muscles, in particular, were verified by an experienced orthopaedic surgeon. Since a single bone segment was defined as the sum of its corresponding soft tissue and bone structures, soft tissue masses could be estimated, based on the subject's mass using regression techniques 69 .
Relevant muscles were implemented in the form of unidimensional Hill-type force elements acting along their anatomical attachment sites. To account for the exertion of physiological forces, the functional units of muscles were subdivided into several force elements to account for wide attachment areas 55 . Deflection phenomena from muscles around the bone were considered by deploying segment-fixed via-points and wrapping 70 . A complete list of muscle structures with their Hill-type model parameters is summarized in the Supplementary Table SI1. Passive forces due to muscular or capsular tissue have been neglected within this work, since they contribute less than 10% to intersegmental torques during comparable activities 60 .
Joint markers and muscle attachment sites-including muscle via-points-adjacent to the flexible femur SEM, required the definition of interface nodes in the FE modeling step as they potentially deform the femoral bone and thus, must be considered during modal reduction. As it concerns the joint markers, the subsets of FE-nodes were defined according to the contacting cartilage surfaces of the femoral head and condyles. As it relates to the muscle attachment sites, the FE-nodes located within a radius of approx. 0.5 cm around the respective attachment site were defined as the interface node's subset of respective coupling nodes. The choice of the radius of 0.5 cm results from FE model's element size of 3.8 mm and the fact that the attachment site area is accounted for by the subdivision of the whole muscle into functional subunits as defined by Carbone et al. 55 and not the area of each subunit itself.
Inverse kinematics. The two-leg squat motion scenario as available from the SimTK dataset 45 comprises the trajectories of reflective skin markers that were applied to the male subject's body at prominent anatomical landmarks to track its motion. The subject performed three squat cycles with a self-selected speed. Utilizing the kinematic topology of the NMS model as defined above, the marker motions and the resulting overall motion of the rigid skeletal multibody system were converted into trajectories of the minimal coordinates in joint space β β β .. , , resulting in a total of seven relative joint coordinates of the kinematic chain. The experimentally obtained joint The musculoskeletal multibody model as described in detail in continuous text. Since the dynamic motion is muscle driven, spring-damper force elements restricted the motion of the pelvis to the sagittal plane and thus, resembled the reaction forces of the lower-left extremity for the two-leg squat. The femur superelement was incorporated using a master node for each joint marker, each muscle attachment site and each via-point, respectively. (2019) 9:6928 | https://doi.org/10.1038/s41598-019-43028-6 www.nature.com/scientificreports www.nature.com/scientificreports/ trajectories were then used as the input for the inverse dynamics and computed muscle control to generate coordinated muscle forces and in turn, the desired squat motion.
Inverse dynamics and computed muscle control. In our NfMBS featuring dynamic stress and strain computation, we deployed a computed muscle control (CMC) algorithm 71  is used for input-output linearization of the NMS system. On the right-hand side of (6), the overall inertial torque which comprises the inertia torque β ⋅M( ) appears due to the acceleration β exp and the torque τ ⋅ ( ) cc which are due to Coriolis and centrifugal effects and the torque due to gravity τ ⋅ ( ) g . The vector r is a set of absolute coordinates that describe the actual spatial positions of the bodies, relative to a fixed inertial reference system  0 . On the left-hand side of (6) the torques due to active muscle forces τ ⋅ ( ) m appear. The NMS system dynamics is then superposed with a generic feedback control, K e K e e e with , , the tracking errors will converge to zero in a critically-damped manner.
The resulting net torques τ m around the respective DoF is then distributed over the available musculotendon actuators ⋅ f ( ) to generate muscle forces. Therein, n is the number of muscles, ⋅ f ( ) i lv, is the force-length-velocity relation with muscle length s and muscle contraction velocity  s kinematically expressed in terms of the joint angles β and joint angular velocities β  , A i is the physiological cross-sectional area, σ i is the maximal isometric muscle stress, and u i is the muscle's unit direction vector.
For the sake of simplicity, the force-length-velocity factor was set to ⋅ = f ( ) 1 i lv, as properties of the activated muscle structures and the activation dynamics itself have little impact on the prediction of muscle forces 72 . The muscle force element is then defined as an actuator whose force generating capacity depends on its theoretical maximum force σ A i i , its activation level a i , and primarily on the moment arm. Note that we used the multibody formalism for rigid bodies, as we do not expect the bone deformations to have any impact on the overall system motion or the muscle force computation.