## 1 Introduction

Three dimensional (3D) printing is an emerging technique of layer-by-layer additive manufacturing technology, with growing interest in medical applications (Rengier et al., 2010). This is because 3D-printed phantoms (or prototypes) can provide precise mimicking of patient organs at an acceptable price and time cost (Ventola, 2014). Such phantoms can be extremely helpful for doctors to practice and be proficient in surgical procedures (Chen et al., 2018b) as well as personalized pre-surgical planning (Qian et al., 2017). However, one limitation is that the mechanical property of printed phantoms is completely different from biological tissues (Raghavan et al., 1996). Currently, the state-of-the-art approach is to print sub-structure, or metamaterial (see Figure 1), to mimic the mechanical properties of biological tissue (Wang et al., 2016a, b). However, such a mimicking procedure requires careful tuning, which may take days or even weeks to perform (Chen et al., 2018a). This greatly limits the medical applicability of tissue-mimicking metamaterial, since surgery timing is a critical factor for outcome success. In this paper, we propose a novel kriging model for emulating mechanical response (functional output) over the design space of metamaterial geometry (functional input), which can be used for efficient tissue-mimicking optimization in practical turnaround times.

There are two reasons why a statistical surrogate model is needed for this tissue-mimicking problem. Firstly, exploration of the design region for metamaterial structure using only physical experiments (i.e., by 3D-printing each phantom, and physically testing its mechanical response) is impractical in practice, due to high material costs and lengthy experimental time. Because of this, one can afford to test only a handful of metamaterial designs, which is not enough to learn the complex functional mechanical response. Secondly, while numerical simulations (e.g., finite element analysis, see Zienkiewicz et al., 1977

) provide a reliable estimate of mechanical response, such simulations can again be quite time-intensive to run. Here, the finite element run for a single metamaterial geometry requires at least 30 minutes of computation; because of this, the time required to explore the large (functional) metamaterial space via only simulations would be far beyond the timeframe available for urgent surgical procedures. To this end, a surrogate model (or emulator, see

Santner et al., 2013) provides an effective way to survey the large design space and perform tissue-mimicking optimization in acceptable turnaround times.The proposed emulator utilizes a popular technique called kriging (Matheron, 1963)

, which models the unknown simulation output via a Gaussian process (GP). Kriging is widely used in computer experiment modeling for its interpolating property, and the fact that both the predictor and its uncertainty have closed-form expressions

(Santner et al., 2013). One challenge for the current problem is that the output of interest (mechanical response, i.e., the deformation at different load situations) is functional in nature. To deal with functional outputs, many works in the literature employ some form of reduced-basis modeling, including principal component models (Higdon et al., 2008; Mak et al., 2018) and wavelet models (Bayarri et al., 2007), where only the weights of the bases need to be modeled. Other works utilize a co-kriging framework (Stein and Corsten, 1991) for modeling a discretization of functional output – such an approach typically works well when the output is highly correlated in the functional domain (Banerjee et al., 2014). For computational tractability, separability is typically assumed on the co-kriging correlation structure (Hung et al., 2015). Another challenge is that the inputs (metamaterial geometry) are of functional form as well. Several modeling techniques have been proposed for functional inputs in the functional data analysis literature (see, e.g., Ramsay, 2005), including varying-coefficient models (Fan and Zhang, 2008) and historical functional linear models (Malfait and Ramsay, 2003). However, the literature on kriging with functional inputs is scarce. For the specific case of time-series inputs, Morris (2012) proposed a kriging model with a covariance function depending on time order, and Muehlenstaedt et al. (2017) used B-spline basis functions to reduce functional inputs to weights for modeling. These kriging models, however, do not incorporate spectral information of functional inputs, which is known to be important in many engineering problems, including the tissue-mimicking application at hand.We propose in this work a new kriging emulator which makes use of the so-called spectral-distance – the Euclidean distance between two functional inputs in spectral domain – within the correlation function. This spectral-distance is motivated by the physics of the tissue-mimicking problem. On one hand, for two input geometries which are the same except for a translation shift, the mechanical properties for such geometries are the same. Our spectral-distance (or SpeD) correlation captures this physics by assigning a correlation of 1 (i.e., perfect correlation) for these two inputs. On the other hand, previous studies (Wang et al., 2016a; Chen et al., 2018a) showed that spectral differences in metamaterial geometries contribute greatly to variations in mechanical property. The SpeD correlation directly models this physics by measuring the distance between two functional inputs in the spectral domain. Another novelty of our model is the use of sparse priors on both the input spectral coefficients and the output co-kriging covariance matrix. These priors allow the emulator to learn important physical information (e.g., dominant input frequencies, output curve properties), which can then be used for effective tissue mimicking in urgent surgical procedures.

The paper is structured as follows. Section 2 gives an overview of the tissue-mimicking problem and the underlying computer experiment. Section 3 presents the proposed SpeD emulation model and its sparse prior specification. Section 4 outlines the algorithm for parameter estimation. Section 5 investigates the emulation accuracy, uncertainty quantification, physics extraction and computation time of the SpeD emulator in a real-world tissue-mimicking case study. Section 6 concludes the work.

## 2 Tissue-mimicking and finite element modeling

We first describe the tissue-mimicking problem (or the metamaterial design problem) and explain the physics of this problem. We then introduce the finite element (FE) analysis as a simulation tool, and provide a brief discussion on experimental design for the FE simulation.

### 2.1 Tissue-mimicking problem

As discussed in Section 1, 3D-printing technology can print patient-specific phantoms with precise geometry (Figure 1 (a)), but these printed phantoms have incorrect mechanical properties compared to the true organs (Figure 1 (b)). The mechanical property under consideration is the stress-strain curve, defined as the internal force (i.e., stress) as a function over the deformation (i.e., strain) under standard external load (Malvern, 1969). The stress-strain curve of the biological tissue typically possesses the property of strain-stiffening, which means the curve is concave upward (see solid blue line in Figure 1 (b)), indicating it becomes stiffer as more load is introduced (Raghavan et al., 1996). However, for 3D-printable material, an opposite property of strain-softening is exhibited (see dotted red line in Figure 1 (b)) due to the plastic-slipping effect and energy dissipation (Hill, 1998).

One way to reconcile this discrepancy is to introduce metamaterial structure (i.e., printed enhancement substructure) within the phantoms (Wang et al., 2016a), to achieve the strain-stiffening property of the biological tissues. Figure 1 (c) shows an example of a metamaterial structure with sinusoidal wave fiber. Here, the stiffer metamaterial structure is designed to have a sinusoidal shape, inside the cuboid matrix of soft material. In this work, we treat the shape of the metamaterial fiber (assumed to have uniform diameter) as the functional input for our SpeD model. Our goal is to mimic the target mechanical property of human tissues or organs, by carefully choosing the shape of the metamaterial fibers. Figure 1 (d) shows a printed “tissue-mimicking” aortic valve with the optimal metamaterial structure.

### 2.2 Finite element modeling and experimental design

In this work, FE modeling is used to simulate the output stress-strain curve of a given metamaterial design. FE is a frequently-used method for stress analysis in solid mechanics; it transforms the partial differential equations to their integral form, so that a piece-wise linear formula can be used to approximate the true deformation profile

(Zienkiewicz et al., 1977). The key advantage of the FE simulation, compared to prototype (physical) experiments, is that high accuracy can be achieved with no material cost or human error. In this work, FE simulations are performed using COMSOL Multiphysics. The overall size of the metamaterial cuboid is by by , with physics-based quadratic tetrahedral elements used for meshing. To calculate the stress-strain curve of the metamaterial, one end of the cuboid is fixed while a series of load levels (up to 15% uniaxial deformation) is applied to the other end. The total computation time for one metamaterial design is around 30 minutes on 24 Intel Xeon E5-2650 2.20GHz processing cores. Figures 2 (a) and (b) show, respectively, the stress contour and stress-strain curve obtained by FE simulation.Here, we use a sinusoidal wave structure for designing the training metamaterial geometries, as such a form exhibits the best strain-stiffening property from a recent study (Wang et al., 2016b). The design space has four parameters: the diameter of the metamaterial fiber , and the amplitude , frequency and initial phase of the sinusoidal wave:

(1) |

The experimental design adopted for the sinusoidal coefficients is the maximum projection (MaxPro) design proposed by Joseph et al. (2015). The advantage of the MaxPro design is that space-filling properties are preserved in all possible projections of the design space, which provides good predictions from a GP. A total of metamaterial geometries are used as the training dataset over the design space. The testing design is an 18-run Sobol’ sequence (Sobol’, 1967), also on the sinusoidal coefficient design space. Figure 2 (c) shows the pairwise scatterplots for the training and testing designs. As we will show in Sections 5.1 and 5.2, even with a small training dataset ( samples), the functional stress-strain predictions from our SpeD emulator are quite accurate, and provides noticeable improvements over a standard kriging model with four sinusoidal coefficients as inputs.

## 3 Emulation model

We now present the proposed SpeD emulation model in three parts. First, we introduce the SpeD model for functional inputs, with the simplified setting of scalar outputs. We then extend this framework to model functional outputs using a co-kriging structure, and provide closed-form expressions for prediction and uncertainty quantification (UQ). Finally, we specify the sparse priors for model parameters.

One key theme of this section is the incorporation of known physics within the emulation model. This not only enables good predictive performance in the complex setting of functional inputs and outputs (see Section 5.1 and 5.2), but also provides a data-driven way to learn important physics for tissue-mimicking (see Section 5.3 and 5.4). Table 1 summarizes these known physical properties and the corresponding emulator assumptions, which we discuss throughout the section.

### 3.1 Spectral-distance kriging model

We introduce first the proposed SpeD model for functional inputs , where is the space of function inputs (to be defined later). For simplicity, assume first the case of scalar outputs (functional outputs are introduced later). For the map from functional inputs to scalar outputs, we propose the following SpeD GP model:

(2) |

where is the scalar process mean and

is the process variance. Here,

is the proposed SpeD correlation function, defined as:(3) |

Here, is a distance function to be defined later, is the modulus of a complex number (where is the unit imaginary number), and

is the Fourier transform from the input space of integrable functions,

, to its spectral space . In this paper, we use the following definition of a Fourier transform for input function :(4) |

Physical Properties | Model Assumptions |
---|---|

Translation-invariance property for inputs | |

(Chen et al., 2018a) | Spectral-distance correlation function |

Frequency information is important | |

(Chen et al., 2018a) | Frequency domain weight function |

Only a few important frequencies | |

(Chen et al., 2018a) | Sparsity in frequency coefficients |

Separability of geometry and mechanics | |

(Fung, 2013) | Separable co-kriging structure |

Known structure of mechanical response | |

(Callister and Rethwisch, 2011) | Sparsity in co-kriging covariance matrix |

Similar to the scale-parametrized distance function in the Gaussian correlation (which is widely used for GP emulation of computer experiments, see Santner et al., 2013), we will use the following scale-parametrized distance function in the spectral domain:

(5) |

Here, is a weight function in spectral space, with a larger value of indicating greater importance of frequency in SpeD correlation. In contrast to the standard Gaussian correlation, we assign importance to each frequency component of a functional input, rather than to each input variable. Plugging (5) into (3), our SpeD correlation can be written as:

(6) |

We will show later in Theorem 1 that is indeed a positive semi-definite kernel.

One advantage of the SpeD correlation is that it can directly capture known properties of the tissue-mimicking problem. First, recall that the translation-shifting property of Fourier transform: for any , if , then

(7) |

For two metamaterial geometries with a shift, i.e., and , we can then show that their outputs are perfectly correlated, i.e.:

(8) |

This is what we call the translation-invariance property of the SpeD correlation. As illustrated in Figure 3, this is a desirable property, since we know from physical knowledge that any translation of the metamaterial geometry does not affect the output mechanical response. The existing functional input models in Section 1, however, do not have this property. Second, it is shown (Wang et al., 2016a, b; Chen et al., 2018a) that the stress-strain curve depends largely on frequency and amplitude but not initial phase in the sinusoidal parametrization (1). One can therefore expect that (i) the Fourier frequencies are significant, and (ii) variations in mechanical response are largely due to differences in frequency intensities . Our SpeD correlation (6) nicely captures both of these properties.

For our tissue-mimicking problem, the specific choice of the Fourier transform with distance of the modulus gives an intuitive parametrization of known physical properties. For other applications, the SpeD correlation (6) can also be used with other spectral transforms (e.g., wavelet transforms) and other distance metrics (e.g., distance). The choice of spectral transform and distance should be made on a case-by-case basis, motivated by prior information from the problem at hand.

The following theorem ensures that the SpeD correlation (6) is a valid positive semi-definite kernel.

###### Theorem 1.

The SpeD correlation function in (6), is a positive semi-definite kernel, i.e.:

(9) |

holds for any and any distinct functions .

The proof is provided in Appendix A. This positive semi-definite property guarantees the validity of as a proper correlation function to use for GP modeling. Note that is not (strictly) positive-definite, in that an equality in (9) does not imply for all . This can be seen by setting all input functions to be the same modulo a translation shift; the resulting correlation matrix then becomes a matrix of ones, which is clearly not positive definite. The fact that is not positive-definite is not an issue here. This is because by the projection property of the MaxPro design (Joseph et al., 2015), all training input functions are distinct even after translation shifts.

### 3.2 Spectral-distance co-kriging model

For the tissue-mimicking problem, the output – the stress-strain curve – is of functional form as well. Below, we generalize the scalar model in Section 3.1 to account for functional outputs. Denote the functional input as and functional output as , where is the output stress at strain level . For our training dataset of simulated geometries, the functional outputs are discretized into

levels, yielding output vectors

, . We assume the following SpeD co-kriging model on :(10) |

where is the process mean vector and is the corresponding covariance matrix function,

(11) |

Here, is the SpeD correlation kernel in (6), and , is a symmetric, positive definite co-kriging covariance matrix quantifying correlations between different output levels.

Equation (11) implicitly assumes separability in the co-kriging covariance structure. Here, separability means the covariance between output levels observed at different functional inputs can be decomposed as the product of the covariance between output levels and the covariance between functional inputs. This separability assumption is used extensively in the literature for reducing computational complexity (Banerjee et al., 2014). However, such an assumption can be justified by the physics of the tissue-mimicking problem as well. Here, the stress-strain curve output is an inherent property pertaining to the metamaterial structure itself, and the stress value at any specific strain level is not directly related to the metamaterial geometry (Diest, 2013). Thus, the covariance between two strain levels and in the stress-strain curve should be similar between different metamaterial geometries. Separability allows us to efficiently train our SpeD emulator for timely tissue-mimicking in real-world surgical problems.

Next, we derive the equations for prediction and UQ. Let

denote the vector of functional outputs of the whole training set. Using the conditional distribution formula of the multivariate normal distribution, the discretized functional response

at a new functional input follows the multivariate normal distribution:(12) |

where is the Kronecker product, denotes 1-vector of elements, and are the process mean vector and co-kriging covariance matrix, and . After algebraic manipulations, the posterior mean and posterior variance can be written in a more concise form:

(13) | ||||

(14) |

where denotes an identity matrix. Equation (13) can be used to predict (or emulate) the stress-strain curve for a new metamaterial design, while Equation (14) can be used to construct a confidence band for quantifying the uncertainty of this prediction.

### 3.3 Prior specification

We now specify independent priors for the model parameters . Since little information is known on the mean prior to data, we assign a non-informative flat prior . For the weight function , we assign independent exponential priors for at each frequency , where:

(15) |

Here, is a rate parameter for the exponential priors. Similar to the Bayesian LASSO (Park and Casella, 2008), the idea behind (15) is to encourage sparsity in , by assigning large prior mass to regions where . This sparsity is desired for two reasons. First, this allows us to identify dominant frequencies in metamaterial structure which influence mechanical response. Second, sparsity in greatly speeds up the tissue-mimicking procedure with the proposed SpeD emulator, since it screens out inert frequencies from the optimization model, which then reduces the number of evaluations from the emulator. This speed-up is paramount for efficient tissue-mimicking in urgent surgical applications (see Section 5.4).

For the covariance matrix , we assign the following prior:

(16) |

Here, is a rate parameter, and is the element-wise norm. This can be viewed as a sparse prior on (Wang, 2012), corresponding to the widely-used graphical LASSO (Friedman et al., 2008) method for sparse covariance estimation. For our problem, this sparsity can be used to identify important physical couplings in the stress-strain relationship, which can then be used to segment the mechanical response curve in a physically meaningful way (more on this in Section 5.3.1).

## 4 Parameter estimation

In implementation, the functional inputs are also discretized to levels. Let and denote the discretized input vectors for both and , respectively; our SpeD kernel in (6) can be written as:

(17) |

where is the discretized weight vector, and is the -th entry of the discrete Fourier transform for the discretized input . Note that the is symmetric because is real-valued (Sorensen et al., 1987); hence, only the first half of is used in (17).

With this input discretization, we adopt a maximum a posteriori (MAP) approach for estimating the parameters . From the GP model in (10) and (12), the MAP estimation of boils down to minimizing the following negative log-posterior (Santner et al., 2013):

(18) |

Here, is the correlation matrix in (12) with scale parameters , and and are the rate parameters for the sparsity priors in Section 3.3.

From a regularization perspective, the two prior terms and in the negative log-posterior (18) can equivalently be viewed as penalty terms on and , with the rate parameters and corresponding to penalization parameters. In this sense, the parameters and control the degree of sparsity imposed on and , with a larger (or ) resulting in a sparser estimate of (or ), and vice versa. In practice, these penalization parameters can be estimated from the data itself, or specified from the problem at hand. For example, if predictive accuracy of the emulator is the end goal, then and can be estimated based on cross-validation techniques (Friedman et al., 2001). However, if the extraction of important physics is desired, then and can be set so that a desired number of physical features can be learned. We will return to this in Section 5.3.

Consider now the MAP optimization in (18) for fixed and . We will use the following blockwise coordinate descent (BCD) optimization algorithm, described below. First, assign initial values for , and . Next, iterate the following three steps until the convergence is achieved: (i) for fixed parameters , compute the correlation matrix and the process mean using closed-form expressions (see Santner et al., 2013 for details); (ii) for fixed GP parameters and , optimize for using the graphical LASSO algorithm (Friedman et al., 2008); and (iii) for fixed mean and covariance matrix , optimize for using the L-BFGS algorithm (Liu and Nocedal, 1989). The full optimization procedure is provided in Algorithm 1. Since (18) is a non-convex optimization problem, the proposed BCD algorithm only converges to a stationary solution (Tseng, 2001). Because of this, we suggest performing multiple runs of Algorithm 1 with random initializations for each run, then taking the converged estimates for the run with smallest negative log-likelihood. These runs should be performed in parallel if possible, to take advantage of the parallel computing capabilities in many computing systems.

## 5 Emulation results

In this section, we present the numerical performance of the proposed model for tissue-mimicking. This is presented in four parts. First, we compare the prediction performance of the proposed SpeD emulation model with two baseline emulation models. Second, we provide a comparison of the uncertainty quantification from these three emulation models. Third, we analyze the physical properties learned via sparse priors on and . Finally, we demonstrate the usefulness of the fitted model for mimicking human aortic tissue.

### 5.1 Prediction accuracy

As mentioned in Section 2.2, the proposed SpeD model is fitted using the training data of FE simulations. The input function is discretized to parts, which we denote as a vector . In this specific tissue-mimicking problem, the diameter of the metamaterial enhancement (assumed to be uniform over the whole functional curve ) is also important. To account for this extra design variable, we use the following separable correlation function :

(19) |

where is the discretized SpeD kernel in (17) and is the scale parameter for diameter . Meanwhile, the output function is discretized to parts, which we denote as a vector .

From the underlying physics of the stress-strain relationship, we know that (i) both the strain and stress are positive, and (ii) the stress is zero when strain is zero (this is known as the free-standing state, see Malvern, 1969). To account for (i), a standard log-transformation of the functional outputs is performed prior to modeling and parameter estimation, and the final results are transformed back to ensure the predicted stress is always positive. To account for (ii), we expect a boundary condition, to hold (approximately) for the emulation predictions, since all stress-strain curves in the training data follow this constraint.

For comparison, we also fit two different emulators are fitted as baseline methods, using the same dataset. Both methods use the same co-kriging framework for discretized outputs . The inputs of the first emulator are the parameters from the sinusoidal wave design , which represents the diameter of the metamaterial fiber, amplitude, period and initial phase of the sinusoidal wave (see Figure 1 and Equation (1)). This emulator uses a GP model with correlation function:

(20) |

The same correlation function (with scalar output) is used in Chen et al. (2018a). We refer this model as a feature-based co-kriging model. The second emulator also assumes a GP model with correlation function:

(21) |

This correlation (21) is essentially the Gaussian correlation function, with distance taken to be the -distance between input functions. A similar correlation function is used in Morris (2012) for time-series inputs, with additional dependencies on time order. We refer this model as a -distance functional co-kriging model.

#### 5.1.1 Predicting stress-strain curve

To test the performance of the proposed emulator, we compare the predictions of stress-strain curves (using Equation (13)) for the metamaterial designs from the testing set (see Section 2.2). Figure 4 shows the emulated stress-strain curves for two test metamaterial geometries, along with the true stress-strain curve (ground truth) from FE simulations. To quantitatively measure the difference between the predicted and true curves, we use the following mean absolute relative error (MARE) metric:

(22) |

where is the strain level, is the stress at strain from FE simulation (ground truth), and is the predicted stress from the emulators. Table 2 reports the MARE values for the two test geometries, using the three considered emulators. The proposed SpeD emulator appears to perform very well in all of the testing cases, in that it achieves noticeably lower MARE than the two existing emulators. Figure 5 shows the boxplot of MARE over the whole 18 testing cases for the three emulators. For the proposed emulator, the mean MARE over the whole testing set using proposed emulator is , which is noticeably smaller than the baseline emulators with mean MARE of and , respectively. The SpeD emulator therefore performs better than the two baseline emulators, in terms of predicting the true stress-strain output curve. This is not surprising, since our model captures known physical properties of the tissue-mimicking problem (see Table 1).

#### 5.1.2 Predicting physical characteristics

In addition to predicting stress-strain curve at untested metamaterial designs, engineers are also interested in predicting key physical characteristics. An accurate prediction of these physical characteristics can be as important as emulating the stress-strain curve itself, because it provides interpretability to the black-box emulation model. Two important physical characteristics to predict here are (i) the elastic modulus of the stress-strain curve, and (ii) the classification of material type as strain-stiffening or strain-softening. For (i), the modulus, i.e., the slope of the stress-strain curve at different strain levels, can be interpreted as the stiffness or hardness of the material (Raghavan et al., 1996). Here, we are interested in the elastic moduli and at strain levels and , respectively, where

; this allows us to evaluate the elastic moduli prediction over a wide range of strain levels. For (ii), we wish to classify the stress-strain curve as

strain-stiffening or strain-softening; this is particularly important given the goal of mimicking biological tissues (see Section 2.1). One way to classify is to use the curvature of the stress-strain curve, which can be approximated by the slope of the two moduli, . Assuming no fluctuations in the strain region (Malvern, 1969), a positive curvature suggests a strain-stiffening property is present (due to increasing moduli), while a negative suggests a strain-softening property is present. Figure 6 (a) visualizes these two physical characteristics from a stress-strain curve.SpeD | Feature-based | -distance | |
---|---|---|---|

True positive | 12/12=100% | 11/12=91.7% | 7/12=58.3% |

True negative | 6/6=100% | 5/6=83.3% | 5/6=83.3% |

Accuracy | 18/18=100% | 16/18=88.9% | 13/18=72.2% |

We now compare three emulators (SpeD and baselines) for predicting the moduli and material type. The moduli and , calculated from the emulated stress-strain curves, are compared with the moduli and from FE simulation. The boxplots of the absolute relative error, i.e., , are shown in Figures 6 (b) and (c). The proposed SpeD model out-performs both baseline emulators here, in that it has smaller median errors and lower error variability over the test set. For classification, the predicted curvature , computed from the emulated curves, are compared with the true curvature from FE simulation. Table 3 shows the correct classification rates for the three emulators. The SpeD model has a perfect classification accuracy: it identified the correct strain-softening/-stiffening property for all 18 test geometries. On the other hand, the feature-based model and the -distance model achieves only a and classification accuracy rate, respectively. One reason why the proposed SpeD model can better capture these physical characteristics (compared to existing emulators) is because it directly incorporates the underlying physics via the SpeD correlation function.

### 5.2 Uncertainty quantification

Particularly in healthcare applications, the quantification of predictive uncertainty can be as important as the prediction itself. For the proposed model, Equations (13) and (14) can be used to construct 90% pointwise highest posterior density predictive intervals (HPD-PIs) for the emulated stress-strain curves. Figure 7 shows the HPD-PI for the three emulation models. One observation is that there is little predictive uncertainty at zero strain, with uncertainty increasing as strain levels increase. The observation is consistent with the physical intuition in Section 5.1: the stress always equals to zero when strain equals zero, i.e., no force at free-standing condition. The increasing uncertainty for higher strain levels may be due to the log-transformation of the functional output.

Of the three emulators, the feature-based model has the narrow HPD-PIs, which is expected since it assumes a sinusoidal parametric form for the functional inputs. However, this narrow HPD-PI fails to cover the true stress-strain curve (see both middle plots in Figure 7) in many of the test cases! For the other two methods, the HPD-PIs cover the true curve from simulations, with the HPD-PIs for the proposed model noticeably narrower than that for the -distance model. In this sense, the SpeD model provides better uncertainty quantification over the other two models: it gives reliable coverage of the true stress-strain curve, with relatively low predictive uncertainty. The reasons for this may be two-fold: (i) the SpeD correlation directly captures the physics of the tissue-mimicking problem, which can be thought of as an additional source of data, (ii) the sparse priors on spectral coefficients screens out inert frequencies, which also helps reduce predictive uncertainty.

### 5.3 Learning physics via sparsity

In addition to prediction accuracy and uncertainty quantification, the SpeD emulator also provides a data-driven approach to learn important physics. This is achieved via the sparse priors on both the covariance matrix and frequency coefficients .

#### 5.3.1 Segmentation of stress-strain curve

We first analyze the important correlations selected by imposing sparsity on (the inverse of) the co-kriging covariance matrix . Recall from Section 3.3 that the graphical LASSO is used to extract significant correlations in . Setting the penalty parameter such that 20% of the entries are deemed significant, Figure 7 (a) visualizes the selected (important) covariances in . Recall that the -th entry of represents the covariance between variables and , conditional on all other variables. We see that the stress-strain curve can be roughly segmented into three regions: small strain (from to ) with little conditional correlation, medium strain (from to ) with high conditional correlation and large strain (from to ) with moderate conditional correlation.

These three regions nicely corroborate known physical properties in material strength (Malvern, 1969; Martin et al., 1998), where the mechanical response of the soft bio-mimicking material can also be divided to three regions: the linear region, the elastic region and the yield region (see Figure 8 (b)). We first look into the linear region. From the inverse covariance matrix in Figure 8 (a), the graphical LASSO selected no significant conditional correlations. This can be explained as follows: from material strength research, the stress-strain curve is known to be almost linear within the first region. Hence, any two observations can determine the underlying linear function, which explains the lack of significant correlations in within the linear region. For the second region, high conditional correlations are detected by graphical LASSO. This can be explained by the strong strain-stiffening (or strain-softening) property of the mechanical response within the elastic region, which results in steep slopes for the stress-strain curves. For the third region, the graphical LASSO suggests the presence of moderate conditional correlations. Physically, this can be explained by the migration of both strain-stiffening and strain-softening properties from metamaterial straightening.

#### 5.3.2 Learning dominant frequencies

The proposed approach can also learn important frequencies which influence mechanical response, via the sparse priors on the weight function (see Section 3.3). Figure 9 (c) shows the MAP estimate of in the spectral space, where the rate parameter is chosen via cross-validation. We see that shrinks to zero at low and high frequencies, with non-zero estimates only for medium frequencies between to . For these two endpoint frequencies, Figures 9 (a) and (b) show the metamaterial structures with frequencies and , respectively.

The selected frequencies in is also in line with the physical understanding of the problem. For low frequencies (Figure 9 (a)), the fluctuation in metamaterial design is too weak to induce any effect on the stress-strain curve, whereas for high frequencies (Figure 9 (b)), the resulting strong fluctuation in metamaterial leads to mechanical properties similar to a straight fiber (given nonzero diameter , see Chen et al., 2018a). While it is known that different frequencies affect mechanical response in different ways, a strict law is difficult to find for engineers. Here, our SpeD emulator sheds light on the influential frequencies, i.e., from to ; those frequencies should be carefully chosen for metamaterial design. We show next how this identification of important frequencies allows us to greatly speed up optimization for the target task of mimicking human tissue.

### 5.4 Mimicking aortic tissue

We now explore the motivating task of mimicking the mechanical properties of a target tissue with the proposed emulator. Recall that the key challenge in applying 3D-printing techniques to medical application is the discrepancy in mechanical properties between target tissues and printed phantoms. Here, the SpeD model can be used to find a good metamaterial design (both geometry and diameter ) whose stress-strain curve matches the desired mechanical property . This is achieved via the following optimization problem:

(23) |

Here, is the optimal metamaterial geometry, is the optimal fiber diameter, and is the conditional (discretized) stress-strain curve in (12) with diameter and geometry . In words, equation (23) aims to find the optimal metamaterial design whose stress-strain curve from the proposed model (conditional on data) is closest to the target curve in terms of mean-squared error (MSE).

This MSE criterion can be further decomposed as follows:

(24) |

Here, and are the conditional mean and variance of , respectively, and is the trace of the matrix . The first term can be interpreted as trying to minimize the -norm between the emulated stress-strain curve and the target curve. The second term can be viewed as trying to minimize the predictive variance of the emulated curve. Such a decomposition is quite intuitive, since we wish to find a metamaterial design whose emulated curve matches the desired curve, but also has low predictive uncertainty from the emulation model.

One difficulty in solving (23) is that the variable is functional in form, and its discretization , can be too high dimensional to optimize numerically. Here is where the extracted significant frequencies from Section 5.3.2 come into play. Let denote the Fourier coefficients of the eight frequencies which were deemed significant by our model (see Figure 9). Using these coefficients as inputs for optimization (and ignoring the other inert coefficients), we get the following lower-dimensional optimization problem:

(25) |

where is the conditional random vector taking the frequencies as input. While this problem is non-convex, it is much lower-dimensional, and can be effectively optimized using standard quasi-Newton methods (e.g., L-BFGS) and random initializations.

This framework (using our SpeD model) offers significant speedups for tissue-mimicking over the current state-of-the-art methods. To see why, consider first the optimization of (23) using only numerical FE simulations: this requires hundreds of evaluations of the optimization objective function, each of which requires around 30 minutes of computation time. This means tissue-mimicking with only FE simulations can require many days of computation, which is clearly unsuitable for urgent surgical planning (Chen et al., 2018a). To contrast, each evaluation of the proposed emulator requires only seconds of computation, which greatly speeds up the mimicking process. Furthermore, by exploiting sparsity in spectral coefficients, the dimension of the optimization problem reduces from 82 to 9 variables. This dimension reduction greatly cuts down on the number of predictions from the emulator, which yields significant reductions in computation time. Such speed-ups are paramount for performing tissue-mimicking in an accurate and timely manner. Section 5.5 provides a further comparison of timing.

Let us now explore the tissue-mimicking performance of the SpeD emulator in a real-world case study. Figure 10 (a) shows the stress-strain curve of a target aortic tissue (in black) from Hockaday et al. (2012), the stress-strain curve from the proposed mimicking procedure (in red), and the curve from an existing mimicking method (in blue) in Wang et al. (2016b). The latter method performs mimicking using only the four sinusoidal metamaterial features (see Section 2.2). Compared to the existing approach which has an MARE (see Equation 22) of , the proposed emulation approach achieves a smaller MARE of , an error reduction of . This improved tissue-mimicking performance can be seen in Figure 10 (a): the red curve (from the proposed method) closely mimics the desired black curve, whereas the blue curve (from Wang et al., 2016b) overestimates stress at all strain levels. In particular, our method gives much better mimicking in small strain regions – this is important in medical applications due to the relatively small strain deformations in the human body.

There are two explanation for this improved performance. First, the existing mimicking approach in Wang et al. (2016b) is too restrictive, in that it uses only four sinusoidal features and not the full functional form of the input. Second, given a fixed timeframe, the proposed emulation-based approach permits a larger number of objective evaluations via our SpeD model. Figure 10 (b) shows the optimal (discretized) metamaterial design from our emulation-based approach, which is clearly not a sinusoidal function. By considering the broader class of functional inputs as well as allowing for more objective evaluations, our model can identify better metamaterial designs for tissue-mimicking. While there is some discrepancy between the proposed stress-strain curve and the target curve, this may be because of the over-simplified optimization scheme and the limited training data available. Improving this will be the focus of future work.

### 5.5 Computation time

Modeling Step | Computation Time (in minutes) |
---|---|

Parameter estimation | 21.72 |

Prediction & UQ at one input setting | 0.01 |

Tissue-mimicking of a target material | 2.69 |

Another appeal of the SpeD emulator is its computational efficiency. Table 4 summarizes the computation time required for each step of the emulation process, with timing performed on a parallelized system of 24 Intel Xeon E5-2650 2.20GHz processing cores. The computation time required for parameter estimation (with cross-validation on ) is minutes, and once we fit the model, we can predict for multiple settings very quickly ( minutes per setting). To contrast, FE simulations requires 30 minutes for each geometry setting. Because of this, our SpeD emulator can effectively perform the tissue-mimicking procedure using only 3 minutes of computation; this greatly improves upon the standard tissue-mimicking approach with only FE simulations, which may require hours or even days to perform (Wang et al., 2016b) with much poorer mimicking performance (see Figure 10)!

Given the positive results reported in this section, the SpeD emulator can serve as an important tool for furthering the use of 3D printed organ phantoms in medical applications. Our model can provide effective and personalized pre-surgical practicing and planning (Qian et al., 2017; Chen et al., 2019) with dramatically lower costs, which then mitigates risk in complex surgical procedures.

## 6 Conclusion

In this paper, a novel emulation model is proposed for mimicking the mechanical properties of biological tissues. The key challenges in this tissue-mimicking problem are the functional form of both inputs (metamaterial geometry) and outputs (mechanical response). To deal with function inputs, we proposed a new spectral-distance (SpeD) correlation function which directly models the effect of each metamaterial frequency on mechanical response. One appealing feature of the SpeD correlation is the so-called translation-invariance property – this accounts for the fact that two metamaterial geometries, which are equivalent modulo a translation shift, have the same mechanical properties. We then presented a separable co-kriging framework for modeling function outputs. Next, under a sparse prior specification, we proposed an efficient algorithm for MAP estimation of parameters. Finally, we applied the proposed SpeD emulator to a real-world tissue-mimicking study, and demonstrated its effectiveness over existing emulators in (i) emulating and quantifying uncertainty on mechanical response, (ii) extracting meaningful physical insights, and (iii) providing efficient and accurate mimicking performance for human aortic tissue. With the development of multi-material 3D-printing technology, the proposed SpeD emulator can play an important role in furthering the impact of 3D-printing in important biomedical applications in surgery planning and healthcare.

## References

- Banerjee et al. (2014) Banerjee, S., Carlin, B. P., and Gelfand, A. E. (2014). Hierarchical Modeling and Analysis for Spatial Data. CRC Press.
- Bayarri et al. (2007) Bayarri, M., Berger, J., Cafeo, J., Garcia-Donato, G., Liu, F., Palomo, J., Parthasarathy, R., Paulo, R., Sacks, J., and Walsh, D. (2007). Computer model validation with functional output. The Annals of Statistics, 35(5):1874–1906.
- Callister and Rethwisch (2011) Callister, W. D. and Rethwisch, D. G. (2011). Materials Science and Engineering, volume 5. John Wiley & Sons NY.
- Chen et al. (2018a) Chen, J., Wang, K., Zhang, C., and Wang, B. (2018a). An efficient statistical approach to design 3D-printed metamaterials for mimicking mechanical properties of soft biological tissues. Additive Manufacturing, 24:341–352.
- Chen et al. (2018b) Chen, J., Xie, Y., Wang, K., Wang, Z. H., Lahoti, G., Zhang, C., Vannan, M. A., Wang, B., and Qian, Z. (2018b). Generative invertible networks (GIN): Pathophysiology-interpretable feature mapping and virtual patient generation. In International Conference on Medical Image Computing and Computer-Assisted Intervention, pages 537–545. Springer.
- Chen et al. (2019) Chen, J., Xie, Y., Wang, K., Zhang, C., Vannan, M. A., Wang, B., and Qian, Z. (2019). Active image synthesis for efficient labeling. arXiv preprint arXiv:1902.01522.
- Diest (2013) Diest, K. (2013). Numerical Methods for Metamaterial Design. Springer.
- Fan and Zhang (2008) Fan, J. and Zhang, W. (2008). Statistical methods with varying coefficient models. Statistics and its Interface, 1(1):179.
- Friedman et al. (2001) Friedman, J., Hastie, T., and Tibshirani, R. (2001). The Elements of Statistical Learning, volume 1. Springer Series in Statistics, NY.
- Friedman et al. (2008) Friedman, J., Hastie, T., and Tibshirani, R. (2008). Sparse inverse covariance estimation with the graphical lasso. Biostatistics, 9(3):432–441.
- Fung (2013) Fung, Y. C. (2013). Biomechanics: Mechanical Properties of Living Tissues. Springer Science & Business Media.
- Higdon et al. (2008) Higdon, D., Gattiker, J., Williams, B., and Rightley, M. (2008). Computer model calibration using high-dimensional output. Journal of the American Statistical Association, 103(482):570–583.
- Hill (1998) Hill, R. (1998). The Mathematical Theory of Plasticity, volume 11. Oxford University Press.
- Hockaday et al. (2012) Hockaday, L., Kang, K., Colangelo, N., Cheung, P., Duan, B., Malone, E., Wu, J., Girardi, L., Bonassar, L., Lipson, H., Chu, C., and Butcher, J. (2012). Rapid 3D printing of anatomically accurate and mechanically heterogeneous aortic valve hydrogel scaffolds. Biofabrication, 4(3):035005.
- Hung et al. (2015) Hung, Y., Joseph, V. R., and Melkote, S. N. (2015). Analysis of computer experiments with functional response. Technometrics, 57(1):35–44.
- Joseph et al. (2015) Joseph, V. R., Gul, E., and Ba, S. (2015). Maximum projection designs for computer experiments. Biometrika, 102(2):371–380.
- Liu and Nocedal (1989) Liu, D. C. and Nocedal, J. (1989). On the limited memory BFGS method for large scale optimization. Mathematical Programming, 45(1-3):503–528.
- Mak et al. (2018) Mak, S., Sung, C.-L., Wang, X., Yeh, S. T., Chang, Y. H., Joseph, V. R., Yang, V., and Wu, C. J. (2018). An efficient surrogate model for emulation and physics extraction of large eddy simulations. Journal of the American Statistical Association, pages 1–14.
- Malfait and Ramsay (2003) Malfait, N. and Ramsay, J. O. (2003). The historical functional linear model. Canadian Journal of Statistics, 31(2):115–128.
- Malvern (1969) Malvern, L. E. (1969). Introduction to the Mechanics of a Continuous Medium. Number Monograph.
- Martin et al. (1998) Martin, R. B., Burr, D. B., and Sharkey, N. A. (1998). Skeletal Tissue Mechanics, volume 190. Springer.
- Matheron (1963) Matheron, G. (1963). Principles of geostatistics. Economic Geology, 58(8):1246–1266.
- Morris (2012) Morris, M. D. (2012). Gaussian surrogates for computer models with time-varying inputs and outputs. Technometrics, 54(1):42–50.
- Muehlenstaedt et al. (2017) Muehlenstaedt, T., Fruth, J., and Roustant, O. (2017). Computer experiments with functional inputs and scalar outputs by a norm-based approach. Statistics and Computing, 27(4):1083–1097.
- Park and Casella (2008) Park, T. and Casella, G. (2008). The Bayesian lasso. Journal of the American Statistical Association, 103(482):681–686.
- Qian et al. (2017) Qian, Z., Wang, K., Liu, S., Zhou, X., Rajagopal, V., Meduri, C., Kauten, J. R., Chang, Y. H., Wu, C., Zhang, C., Wang, B., and Vannan, M. A. (2017). Quantitative prediction of paravalvular leak in transcatheter aortic valve replacement based on tissue-mimicking 3d printing. JACC: Cardiovascular Imaging, 10(7):719–731.
- Raghavan et al. (1996) Raghavan, M. L., Webster, M. W., and Vorp, D. A. (1996). Ex vivo biomechanical behavior of abdominal aortic aneurysm: assessment using a new mathematical model. Annals of Biomedical Engineering, 24(5):573–582.
- Ramsay (2005) Ramsay, J. (2005). Functional Data Analysis. Springer Series in Statistics, NY.
- Rengier et al. (2010) Rengier, F., Mehndiratta, A., Von Tengg-Kobligk, H., Zechmann, C. M., Unterhinninghofen, R., Kauczor, H.-U., and Giesel, F. L. (2010). 3D printing based on imaging data: review of medical applications. International Journal of Computer Assisted Radiology and Surgery, 5(4):335–341.
- Santner et al. (2013) Santner, T. J., Williams, B. J., and Notz, W. I. (2013). The Design and Analysis of Computer Experiments. Springer Science & Business Media.
- Sobol’ (1967) Sobol’, I. M. (1967). On the distribution of points in a cube and the approximate evaluation of integrals. Zhurnal Vychislitel’noi Matematiki i Matematicheskoi Fiziki, 7(4):784–802.
- Sorensen et al. (1987) Sorensen, H. V., Jones, D., Heideman, M., and Burrus, C. (1987). Real-valued fast Fourier transform algorithms. IEEE Transactions on Acoustics, Speech, and Signal Processing, 35(6):849–863.
- Stein and Corsten (1991) Stein, A. and Corsten, L. (1991). Universal kriging and cokriging as a regression procedure. Biometrics, 47(2):575–587.
- Tseng (2001) Tseng, P. (2001). Convergence of a block coordinate descent method for nondifferentiable minimization. Journal of Optimization Theory and Applications, 109(3):475–494.
- Ventola (2014) Ventola, C. L. (2014). Medical applications for 3d printing: current and projected uses. Pharmacy and Therapeutics, 39(10):704.
- Wang (2012) Wang, H. (2012). Bayesian graphical lasso models and efficient posterior computation. Bayesian Analysis, 7(4):867–886.
- Wang et al. (2016a) Wang, K., Wu, C., Qian, Z., Zhang, C., Wang, B., and Vannan, M. A. (2016a). Dual-material 3D printed metamaterials with tunable mechanical properties for patient-specific tissue-mimicking phantoms. Additive Manufacturing, 12:31–37.
- Wang et al. (2016b) Wang, K., Zhao, Y., Chang, Y.-H., Qian, Z., Zhang, C., Wang, B., Vannan, M. A., and Wang, M.-J. (2016b). Controlling the mechanical behavior of dual-material 3D printed meta-materials for patient-specific tissue-mimicking phantoms. Materials & Design, 90:704–712.
- Zienkiewicz et al. (1977) Zienkiewicz, O. C., Taylor, R. L., Zienkiewicz, O. C., and Taylor, R. L. (1977). The Finite Element Method, volume 36. McGraw-Hill London.

## Appendix A Proof of Theorem 1

For any and function , we have

(26) | ||||

(27) |

Note that the Fourier transform , where is the space of integrable functions , has a unique inverse . Denote the standard Gaussian kernel as ,

(28) |

where is the space of integrable functions . Since is a positive definite kernel, for the selected and in Equation (26), and any function , we have,

(29) |

Now let , where . This is possible because and therefore . Thus, we have

(30) |

In other words,

(31) |

i.e., the proposed SpeD correlation function is positive semi-definite.