Characterization of Spatially Graded Biomechanical Scaffolds

Over the last two decades, there has been an increasing interest in using bio-engineered constructs for tissue engineering, for everything from making artificial ligaments to enhancing the body's ability to heal bone, skin, and muscle [1–6]. Nanofibrous scaffolds can be used to transform this research; by carefully constructing these scaffolds, model systems can be made strong and sufficiently pliable to mimic load directions placed on bone and connective tissue, can degrade naturally in the body, and can be engineered to be porous and hydrophilic to promote cellular response and allow the introduction of growth factors through the scaffolds themselves [7–9].

Characterizing these scaffolds is not a trivial proposition; on a microstructural level, there are numerous qualities of importance that can affect the mechanical response, including fiber size, topography, and alignment, as well as the size and orientation of pores. Characterization at a larger scale provides the mechanical properties of the scaffold as a whole, which are critical for ensuring that it is able to bear the mechanical stresses inherent in a given application [6]. There are a variety of methods, such as scanning electron microscopy (SEM) and mercury porosimetry, for assessing morphology and assorted structural properties; but, for characterizing mechanical properties, the field is largely dominated by tensile testing, with compression testing being employed for special cases and some in vitro studies [6,10,11]. Recently, a number of studies have sought to characterize the mechanical properties of individual electrospun fibers, using either atomic force microscopy (e.g., Refs. [12] and [13]) or single-fiber tensile testing [14–16]. Additionally, complex constitutive models have been used to directly link scaffold microstructure and overall mechanical properties, but these tend to be tailored to highly specific applications and materials (for examples, see Refs. [12] and [17]).

Despite these advances in the microstructure-informed scaffold mechanical behavior, the bulk of the work reported in the literature uses scaffolds with properties that are designed to be uniform throughout, typically consisting of a single layer of nonwoven fibers [17–21]. While these scaffolds can provide a useful environment for the growth of new cells and are able to mimic simple structures, they are poorly suited for materials with graded properties, such as those encountered at interfaces between dissimilar tissue, including the tendon-to-bone transition, the myotendinous junction, and ligament-bone entheses. These structures are inherently nonuniform, with significant gradations in qualities like mineralization, fiber alignment, and modulus throughout [22]. Successfully mimicking these structures requires not only more advanced fabrication techniques to replicate these material property gradients, but also appropriate characterization techniques that are able to discern spatial variations in material properties—which can be very challenging [20,23–26]. While some methods for morphologic characterization, like SEM, still work for spatially graded scaffolds, standard tensile characterization methods are wholly deficient for spatially varying mechanical properties like modulus, as these properties are homogenized throughout the specimen geometry. Thus, a tensile characterization technique, able to determine material properties and how they vary throughout the specimen geometry, is highly desirable, and would be of great interest to bio-engineering applications at disparate tissue interfaces.

In 2009, Li et al. created a scaffold with graded properties, characterized its microstructure, and measured how its mechanical properties were distributed throughout [27]. In their study, the tendon-bone insertion site was mimicked by creating an electrospun scaffold with a continuous gradient in mineralization, and a variety of methods, including SEM and X-ray spectroscopy, were used to examine its microstructure. They then applied a tissue-staining agent to create a speckle pattern on the sample, and performed tensile tests in conjunction with image correlation techniques to record displacements at numerous points on the sample surface and generate images of the local strain throughout the sample as it underwent elongation. Their results showed a clear correlation between the gradients in mineralization and modulus, as expected. This type of strain imaging gives a sense of the distribution of the modulus throughout the scaffold, greatly enhancing the usual tensile testing procedure. However, this approach cannot appreciate the nonuniformity of stress throughout heterogeneous tissues, and is particularly sensitive to noise due to its numerical differentiation approach to extract strain.

Herein, we aim to address these limitations and extend our characterization capabilities through elastography, by creating our own set of graded nanofibrous scaffolds and solving for the distribution of their mechanical properties, namely, the Young's modulus, by solving an inverse elasticity problem using displacement data from tensile tests. By reconstructing the modulus directly from the displacement, we omit conversion to strain and ensure that the results are consistent with the laws of elasticity and a given constitutive model [28]. Consequently, though we restrict ourselves to modulus here, this general method also has the advantage of being extensible to other quantities of interest: with a sufficiently detailed constitutive model, we could use these techniques to determine a variety of microstructural properties from displacement fields [29,30]. A number of studies have used ultrasound or magnetic resonance imaging (both in vitro and in vivo) to noninvasively evaluate scaffold properties as they degrade [31–34]. However, to our knowledge, there are no studies that use elastography to bridge the gap between fabrication and mechanical characterization of graded scaffolds, and allow for the assessment of spatially varying material properties.

Multilayered nanofiber mats were created using previously described electrospinning methods [35]. Briefly, a polycaprolactone (PCL) solution, consisting of poly(ϵ-caprolactone) (80,000 MW, BrightChina Industrial, Shenzhen, Guangdong, China) in a 3:1 chloroform-to-methanol solvent (VWR International, LLC, Radnor, PA), was continuously expelled from a glass syringe through an 18-gage needle, using a programmable syringe pump (BS-8000, Braintree Scientific, Braintree, MA). The polymer solution was expelled at a rate of 12 μL/min, through an applied 18 keV field, and collected on a rotating mandrel (15 × 1.5 in.) in a custom electrospinning apparatus. Mandrel rotation speed was controlled using a commercial overhead high-speed stirrer (BDC6015, Caframo Ontario, CA), and maintained at 2000 rpm during deposition to achieve aligned nanofibers. This alignment was characterized by a highly nonuniform distribution of fiber orientation angles across the sample, with most fibers in the neighborhood of a particular orientation, similar to that seen in our previous electrospun aligned PCL nanofibrous mats [35]. To create scaffolds with spatially graded properties, we utilized a seven-layer design (Fig. 1) wherein sections of select layers were masked to prevent nanofiber accumulation, using card stock strips positioned above the mandrel during deposition (Fig. 2). Because these masking strips were not in contact with the mandrel surface, and thus, did not contact the specimens directly, their removal did not affect the material. In the present scaffold design, continuous layers were alternated with interrupted (masked) layers such that the resulting multilayer structure was symmetric in specimen thickness, with continuous (unmasked) layers as the initial base, and final layers. This design produced scaffolds with uniform continuous surfaces and aligned nanofiber topography, with 80% of fiber diameters between 200 and 425 nm (Fig. 1).

Following fabrication, the multilayered nanofiber mats were removed from the mandrel, and cut into individual material test specimens using a 3.6 × 0.6 cm rectangular punch (C.S. Osborne, Harrison, NJ). Specimens were obtained such that fibers were aligned to the long axis of the sample, and the masked regions were centered within the specimen gage length (i.e., samples were most compliant in the center, and stiffest at the ends) (Fig. 3). Each individual test specimen was then mounted in an oak-tag I-frame sample holder (Figs. 4(a) and 4(b)), and an anisotropic, high-contrast speckle pattern was applied to enable previously described digital image correlation (DIC) measurement [35]. Due to the high degree of fiber alignment, we expected these structures to display anisotropic behavior during tensile testing.

In addition to preparing samples for tensile testing, during one of the electrodeposition sessions, we created six additional samples for destructive thickness measurements. The thickness samples (n = 6) were analyzed using a metallurgical microscope, which allowed void spaces to be removed without deforming the nanofiber layers. Macroscopically, these samples displayed a fairly uniform macroscopic thickness, comparable to the thickness of the seven-layer region. Metallurgical microscopy yielded thicknesses of 60.93±2.55 μm, 21.35±1.07 μm, and 14.73±0.96 μm for the ends of the samples, the five-layer, and four-layer regions, respectively. Though these thicknesses do not correspond directly to those of the samples used for mechanical characterization, due to the small standard deviations in their thickness measures for each region (thin, medium, thick), they provide an accurate estimation of the graded thicknesses for samples fabricated on that day.

Specimens (n = 4) were mechanically characterized via uniaxial tensile testing with a test resources universal materials test machine, equipped with a 100 N tensile load cell (0.5 N resolution) (Model 500 LE2-1, Test Resources, Inc., Shakopee, MN). In our established method [35], once secured in the pneumatic grips (Fig. 4(c)), the I-beam holder is transected, leaving only the scaffold sample spanning between the grips (Fig. 4(d)), and the specimen is stretched to remove any slack, to a slight positive preload of 0.1 N, corresponding to less than 0.5% strain, which accounts for less than 2% of the maximum applied load. Upon preload, the specimen is stretched to failure at a constant rate (0.2 mm/s), with load, displacement, and time recorded at a rate of 50 Hz. Testing was performed such that the axis of expansion was along the specimen's long axis, referred to hereafter as the axial direction, and the perpendicular axis aligned with the sample width is termed the lateral direction. A two-camera digital image correlation (DIC) system (Correlated Solutions, Inc., Columbia, SC) was employed to measure strains in the material, in a noncontacting manner, using a texture-mapping algorithm. The system was calibrated in multiple planes within the testing volume using a rectangular calibration grid (9 × 9-25 mm). Video was collected at 7 frames/second during tensile testing, and analyzed with VIC-3D 2010 dic software (Correlated Solutions), to produce a texture map, treating the solution as a continuum of displacement data across the sample. A virtual extensometer was created along the centerline of the sample length (axial direction) using vic2d software (Correlated Solutions) to determine the axial engineering strain of the sample within its gage length [35]. In order to make the data compatible with our inverse elasticity problem solver, which uses the finite element method, we discretized the specimen geometry using a 15 × 40 rectangular virtual grid, conforming as closely as possible to the boundaries of the calibration region, approximately 5 × 25 mm. The 3D positions for every node of the virtual grid were exported at each time-step to determine the measured displacements. This resulted in a series of incremental frame-to-frame displacements at each point on this new mesh, with the corresponding load cell and grip-to-grip displacement values providing the boundary conditions at each time-step.

Principal component analysis (PCA) allows us to take a large set of measurements and decompose it into a number of orthogonal modes equal to the number of measurements [36,37]. When each measurement is a repetition of the same experiment, e.g., ten different measurements of the length of a sample, or 20 different tests of a component's resistivity, these modes can be used to separate noise from signal. In this case, the lower modes, which correspond to the larger singular values, are considered to be signals, and the high modes are considered noise. Here, we use each individual frame of the displacement data as a different incremental measurement. As long as we restrict ourselves to the small-strain regime, taken to be below 5% strain, the linear behavior of the material will yield approximately the same displacement for each incremental stretch, allowing us to treat each frame as a repetition of the same measurement in our analysis. Our 5% estimate for the cutoff of the linear regime is based on our prior studies in similar electrospun PCL scaffolds [35]. We exclude frames of data collected in the nonlinear, strain-stiffening region, or “toe” region, of the stress–strain curve (Fig. 5), and include only frames for which the average local strain (determined by grid size) across the sample is below 3.5%; this limit is somewhat arbitrary, but justified based on further PCA analysis in Sec. 3.2. PCA is performed on the lateral and axial components of displacement separately, producing modes that represent characteristic motions of the sample in each direction.

Each of these modes generated by PCA accounts for a certain percentage of the total variance present in the data. Looking at the largest contributors can give an estimate of the true signal; modes which contribute a disproportionately small amount to the variance can be considered to be noise. Typically, the variance contribution of each mode is plotted and they are ordered from greatest to least. Of all the modes, we retain the first few dominant ones that together contribute a significant portion of the total variance, around 99%. Thereafter, we project the measured displacement onto these modes to arrive at the denoised estimate of the displacement field.

The first term in π represents the mismatch between a predicted displacement field u and the measured displacements $ũ$⁠. The predicted displacements are constrained to obey a model of linear motion, the details of which are given in Sec. 2.5. T is a matrix used to apply weights to the data, in case the quality of data is dependent on the direction of measurement, as is the case for imaging modalities like ultrasound. As an example, we might account for poor lateral data by using $T=[0.1001]$ to apply a weight ten times greater to data in the axial direction than those in the lateral, thereby forcing the minimization to place more emphasis on the axial data. The R term is a total variation diminishing (TVD) regularization term, depending on the material property distribution μ, as detailed in Refs. [38] and [39]. The overall strength of the regularization is controlled by a regularization parameter α, and E represents the spatial distribution of the Young's modulus. The minimization problem is solved using a quasi-Newton algorithm, which requires calculating the gradient of the objective function with respect to the optimization parameters. This is evaluated efficiently by setting up and solving an adjoint problem [40]. Our approach to solving the inverse problem iteratively proceeds as follows:

1. For a given spatial distribution of material properties, solve a forward elasticity problem to determine the predicted displacements.

2. Compute the objective function, which measures the difference between the predicted and measured displacements. If this is below a given tolerance, stop. Else continue.

3. Solve an adjoint problem that is driven by the difference between the predicted and measured displacement fields.

4. Use the forward and the adjoint problems to determine the gradient of the objective function.

5. Provide this gradient to the quasi-Newton algorithm to determine the updated spatial distribution of mechanical parameters.

6. Go to step 1.

In order to solve the inverse problem, we first consider the properties we expect from the graded PCL specimens. Past studies have shown that similar constructs exhibit anisotropic behavior, due in part to fiber realignment [35,41]. We anticipate behavior to be linear at small strains, but, if realignment is substantial, we may need nonlinear modeling to capture the rotational character of that motion. We also understand that, since the specimens are very thin sheets, we can assume plane stress conditions and work in two dimensions. Taking all this into consideration, to compute multiple components of the Young's moduli, we require a two-dimensional model that can account for anisotropy and accommodate variable axes of anisotropy, and is dependent on Poisson's ratio. To use this model, we need an estimate for Poisson's ratio as well as multiple two-dimensional displacement measurements from experiments stretching the sample along different axes. However, performing multi-axial loading without irreversibly altering the scaffold microstructure from test to test is very difficult, and our previous work has shown that the complex behavior of these structures makes Poisson's ratio challenging to assess [35]. Designing experiments to account for these difficulties and perform suitable multi-axial loading is beyond the scope of this study, so we limit ourselves to a simpler experiment and model.

Our experiment, as described in Sec. 2.2, consists of a single tensile test, stretching samples uniaxially in the direction of fiber alignment. We note that in this case, only the lateral motion is influenced by Poisson's ratio, and the axial motion is governed entirely by the corresponding axial component of the Young's modulus, E₁₁. Thus, we can reduce the two-dimensional model to a one-dimensional model, wherein isotropy holds and we recognize that the only material parameter for which we can solve is E₁₁. In spite of this reduced dimensionality, the inverse problem we have described is driven by the solution of a forward problem in two dimensions; we enforce consistency with our simplified case by letting the lateral component of weighting matrix T be zero. This allows us to choose an arbitrary value for Poisson's ratio, as it will only affect the lateral data, which is ignored by the data-matching term in our optimization.

We note that though this model requires only the axial data, we perform PCA and denoising for both components of displacement in order to interrogate our assumptions about the material behavior and build a framework to do this analysis with more complex experiments and models in the future.

Once the material model and weighting matrix are established, regularization must be addressed. The inverse problem algorithm uses a TVD functional to smooth our reconstructions and ensure that no over-fitting of data occurs. This TVD regularization penalizes large spatial variations in the material parameters without regard to their slope. This makes it an ideal choice for this application, as we expect sharp changes in material properties over small distances due to the spatial gradations engineered into the sample [42]. The strength of the regularization is controlled by a parameter, α, whose optimal value is commonly determined using an L-curve. In general, an L-curve plots a predicted solution against the difference between the predicted and true solutions, for a range of regularization parameter (α) values [43]. We instead use a plot directly comparing the difference between the predicted and true solutions to α itself, referred to by Paynter et al. as an S curve [44]. We look for a point on the bend of the S curve, where the tail of the S begins; this represents a nearly optimal compromise between minimizing the mismatch between measured and predicted results, and avoiding over-fitting the data. We find that an α value of 1 × 10⁻³ is appropriate for all of our samples.

In the course of solving the full inverse problem, we need to solve a forward elasticity problem. In order to produce accurate solutions and ensure that the inverse problem is well-posed, we must prescribe appropriate boundary conditions [45]. In prescribing these, we note that our samples are fixed at both ends in the grips of the tensile test apparatus; one of the grips remains stationary while the other moves, and nothing comes into contact with the sides of the samples. So, we assert that the sides are traction-free and then use displacement boundary conditions on the top and bottom sample edges to match the measured displacement exactly. On the bottom edge, we constrain both the lateral and axial movement with displacement boundary conditions, but on the top edge we constrain only the axial movement. We do this because the bottom edge of the rectangular grid, on which we perform DIC and PCA, is set close to the fixed grip, and we believe it to be strongly constrained. However, the top edge is more free to displace (narrow) laterally, and thus should be less constrained.

Spatially graded samples displayed a characteristic bilinear mechanical response in their stress–strain curves (Fig. 5). This behavior, along with the trends in PCA modes described in Sec. 3.2, allows us to confirm that we are analyzing displacements in the initial linear region for the scaffolds.

For each of the four samples, we identify a series of successive frames of incremental displacement data corresponding to an average cumulative strain of 3.5%; on average, 17 frames are used for each sample. Using these data, we compute the PCA modes for all sets of lateral and axial data. The contribution from each of these modes to the overall variance for a typical sample is shown in Fig. 6, showing that the first two modes in the lateral case, and the first mode in the axial case, contribute significantly to the variance. The total contribution from all the other modes is less than 2.5% in both cases. The dominance of two modes in the lateral data can be interpreted as an indicator of nonlinear behavior, suggesting more complex behavior than in the axial case.

In Fig. 7, we have plotted the first three lateral and axial modes of motion for one particular case. The first shear mode can be interpreted to reflect the realignment of fibers throughout the structure. As discussed in Sec. 2.1, we anticipate that for each sample, the fibers have a dominant direction of fiber orientation. Further, though we have tried to ensure that this direction is aligned with the axis of stretch in the experiment (i.e., the axial direction), the distribution of fibers will create some slight misalignment between these directions. Consequently, during the stretch, the fibers will tend to reorient themselves along the axis of stretch, causing an overall rotation. This can be observed in the first lateral mode in Fig. 7, which corresponds to a net counterclockwise rotation, as we surmised in our modeling considerations in Sec. 2.5.

The second lateral mode corresponds to lateral narrowing, and the third mode is comprised of random motions, as are all higher modes. We note that the first axial mode corresponds to stretch along the axial direction. The second axial mode represents some form of rotation; however, from Fig. 7, we observe that its contribution is three orders of magnitude smaller than that of the first mode. The third axial mode is comprised primarily of noise. These modes are consistent with the early stages of a uniaxial tension test, in which fibers are realigning and the largest displacements are axial.

We would like to emphasize the fact that the axial movement in the first mode is of much higher magnitude than any of the lateral motions. Additionally, only one mode is significant for describing the axial motions, whereas the lateral motion requires two modes; this reinforces our notion that the lateral motion is nonlinear and more complex than the axial motion. To better explore the behavior of the PCA modes and their significance, we plot the contribution to the variance of the first three modes in both the lateral and axial cases, covering a wide range of cumulative strains (Fig. 8).

Turning our attention first to the lateral data in Fig. 8, we see it always has two dominant modes, even at very small strains. We ascribe this to (i) the rotation of re-orienting fibers, which is manifested in the first lateral mode as a shear and, (ii) to the overall extension of the sample. Accounting for both of these is critical to describing the lateral motion in future studies with more complex models. After a short initial period, the second mode, representing lateral narrowing, comes to represent an increasing proportion of the variance. The axial data, on the other hand, tells a simpler story: we see a single dominant mode throughout, indicative of linear behavior. Taking both sets of data into account, we conclude that selecting a strain limit of 3.5% will yield linear axial motion and nonlinear lateral motion.

Returning to the PCA data, we determine the denoised displacement data by projecting the measured displacement onto the first two lateral modes and the first axial mode. The resulting displacement fields for the sample are shown in Fig. 9. Here, for both the lateral and axial displacements, we have plotted the original displacement field, the denoised field, and the difference, which represents our estimate of noise. Based on this analysis, we estimate that the magnitude of the noise, which is approximately equal for both lateral and axial data, is about 1% of the magnitude of the axial data and 10% of the magnitude of the lateral data. This is a substantial difference in signal-to-noise ratio, giving us more confidence in the axial data, both before and after denoising. This is a potential issue that will have to be navigated in any studies that wish to reconstruct the complete two-dimensional properties of these scaffolds.

We produce reconstructions of the relative modulus values throughout the samples; the distribution is solved up to a multiplicative factor. The results for our four samples are shown in Fig. 10.

We can clearly see the spatial gradation in modulus with which the sample was constructed, with light-colored regions of low modulus in the middle and darker regions of high modulus on the ends. This trend is observed in all samples. However, we also observe variability in the modulus distribution between samples. This can be ascribed to either the manufacturing process or the error inherent in our approach to measure displacements and then infer the modulus from it. The average value of the contrast in the modulus between the outer regions of the sample and the central region is reported to be 1.5 (Table 1). We also observe some variability in contrast, around 12%, between samples. The contrast is computed as the ratio of the average modulus values at the ends to the moduli in the center of the sample. In order to evaluate the average value of the E in the central region, a region in the middle, away from the sharp transition in the modulus field was chosen. Whereas for the average value in the outer regions, two distinct regions, one at the top and the other at the bottom, both away from the sharp transitions in the modulus, were chosen.

To get a sense of the contrast in modulus that we should expect between these two regions, we turn to two separate estimates. First, we recall our scaffold design (Fig. 1), which shows that the samples are comprised of seven complete layers on the ends and four in the center, with a five-layer transition zone. The figure implies that the soft region should be between 0.5 cm and 1 cm in length (which is reflected accurately in our reconstructions), and gives us the simplest estimate of modulus contrast of 7/4 = 1.75 between the stiffest and softest regions. Next, we consider our previous study [35] in which we used our electrospinning approach to fabricate sets of PCL scaffolds with uniform thickness, for a range of scaffold thicknesses, and evaluated their elastic moduli [35]. This study showed that, counterintuitively, modulus increased with decreasing specimen thickness; indicating that modulus cannot simply be scaled by cross-sectional area to estimate values for our current samples. Using the results for uniform specimens along with the thickness measurements from Sec. 2.1, we can extrapolate values of 189 MPa, 242 MPa, and 251 MPa for the three regions of our samples, from stiffest to softest. Because the scaffolds are relatively uniform in macroscopic thickness, with most of the difference between soft and dense regions being accounted for by different degrees of fiber packing, we adjust the thickness values accordingly, yielding final modulus estimates of 189 MPa, 85 MPa, and 61 MPa for each region, producing an approximate contrast of three between the stiffest and softest sections.

Returning to Table 1, the contrast for our reconstructions is 1.5; a value below both of these estimates. This is a 15% reduction from the first estimate, computed based on the number of layers in the sample, and a 50% reduction from the second estimate, determined by extrapolating modulus values from previous tests of uniform scaffolds. The first deviation can be explained by the lessening of contrast inherent in the regularization process, but the second is too large to be explained by the regularization or factors in the optimization process. We postulate that the notable differences in fiber packing densities between our graded samples and uniform specimens of similar thickness necessarily change the fiber–fiber interactions (e.g., friction, entanglement), which, in turn, is reflected in changes in the mechanical response significant enough to exaggerate the expected contrast in moduli.

By using PCA to process the displacement data in this problem, we create both a set of denoised data and an estimate of the magnitude of the noise itself. We have shown that the former helps us to enhance the robustness of the inverse problem solution and detect nonlinear behavior, and the success of this method in other studies [46] gives us confidence that PCA is a useful part of our workflow. The set of denoised displacements obtained after PCA can also be utilized to create clean images of strain in the material, where spatial derivatives of noisy fields are often known to cause difficulties. The assessment of the magnitude of the noise via PCA allows us to see how strongly the axial and lateral signals are influenced by the noise, and to adjust hyperparameters, such as the values of weighting matrix T, accordingly. However, after accounting for the noise in the data, potential sources of error remain in our boundary conditions and material model.

Error induced by the boundary conditions likely stems from our treatment of the top and bottom edges of the samples; we require that both components of displacement match the measured data exactly on the bottom edge, and that the axial components match on the top edge. The stiffening effects of these conditions and the influence of the noise, on the lateral data in particular, make this direct enforcement potentially problematic. In future studies, we could avoid assumptions about the boundary conditions altogether by employing special coupled formulations [47,48].

Our constitutive model assumes linear material behavior and isotropy, and does not take into account the fibrous nature of the mats. Both linearity and isotropy are reasonable assumptions for our experiment, given that we work with data in the small-strain regime and only consider one component of the displacement data. Because PCA shows us that the motion at low strains is dominated by fiber realignment, it is clear that we will need to update our models to represent that behavior in order to increase our modeling accuracy. Updating the constitutive model will become much more important in future studies that work with finite strain data, or wish to characterize other material properties; we can build on the work of recent studies that have attempted to capture the effects of fiber orientation and interaction [12,17,29].

Though we use a simple model for this case, our reconstructions have a relatively low variability in contrast, as noted in our earlier consideration of Table 1. Based on previous studies, this level of variability is within what we might expect for this type of inverse reconstruction method [40,49,50]. Additionally, some of this variation may be a product of variability in the manufacturing process, though this is difficult to quantify. The actual spatial variation of modulus across each sample shows rapid transitions between regions of low and high modulus for all cases, and uneven edges for these regions in two of the four cases. Without more detailed experimental characterization of each sample we cannot say that features like the uneven edges are real with absolute certainty, but we do know that total variation regularization tends to produce “staircasing” in results, yielding sharp edges in reconstructions [51,52]. Repeating these tests with an H¹ regularization term produces very smooth transitions, and looking at the scaffold design (Fig. 1) suggests that we might expect to see a two-step transition. Thus, we can assume that the single sharp transition is a result of the total variation regularization combining the two transition steps into one steep divide. To consider whether we could have produced a more detailed result, we examine the finite element mesh used, as shown in Fig. 11, we can see that the sharp transitions in modulus occur within a single element. This leads us to conclude that we are producing results that appropriately represent the granularity of this data. To resolve finer details in these reconstructions, we would require higher resolution in the displacement data, on a finer virtual grid.

We were able to successfully create and characterize spatially graded nanofibrous scaffolds. Four samples were fabricated using electrospun PCL fibers; as material was deposited, the center of each sample was periodically masked to create a gradation of material stiffness mimicking that of natural tendon or ligament. Each of these samples was stretched uniaxially to failure, and a two-camera DIC system recorded displacements across the sample faces during these tests. We used PCA to remove noise from this data and generate a set of clean input displacements for an inverse elasticity algorithm, which was used to infer the spatial distribution of the modulus across each sample.

Our reconstruction of the modulus revealed spatial distributions matching our expectations based on the construction of the spatially graded samples. In future work, we can refine our approach through a combination of experimental and computational changes. Experimentally, we can use more complex tools for characterization of the scaffolds before tensile testing to assess more material properties and make possible quantitative, rather than relative, modulus reconstructions. In order to address anisotropy, we can implement an anisotropic model in our solver and perform tensile tests along multiple axes. To further enhance the utility of our method, this anisotropic model can be built to model fiber-level interactions and properties.

We would also like to acknowledge the invaluable support of Taylor Anderson and Kristen Lee in performing experiments and collecting and formatting DIC data.

• National Science Foundation (CAREER Award CBET-0954990 (DTC); Funder ID: 10.13039/100000001).

• Rensselaer Polytechnic Institute's Presidential Scholars Fellowship (Funder ID: 10.13039/100007092).

DIC =

digital image correlation

E =

Young's modulus

PCA =

principal component analysis

PCL =

polycaprolactone

R =

regularization function

T =

weighting matrix

TVD =

total variation diminishing regularization

u =

predicted displacement field

$ũ$ =

measured displacement field

α =

regularization parameter

π =

objective function

Ω =

domain; region of interest on a sample