Hierarchical Bayesian constitutive model selection for high-strain-rate soft material characterization

Abstract

The high-fidelity characterization of soft, tissue-like materials under ultra-high-strain-rate conditions is critical in engineering and medicine. Still, it remains challenging due to limited optical access, sensitivity to initial conditions, and experimental variability. Microcavitation techniques (e.g., laser-induced microcavitation) have emerged as a viable method for determining the mechanical properties of soft materials in the ultra-high-strain-rate regime (higher than 10³ s⁻¹); however, they are limited by measurement noise and uncertainty in parameter estimation. A hierarchical Bayesian model selection method is employed using the Inertial Microcavitation Rheometry (IMR) technique to address these limitations. With this method, the parameter space of different constitutive models is explored to determine the most credible constitutive model that describes laser-induced microcavitation bubble oscillations in soft, viscoelastic, transparent hydrogels. The target data/evidence is computed using a weighted Gaussian likelihood with a hierarchical noise scale β, which enables the quantification of uncertainty in model plausibility. Physically informed priors, including range-invariant, stress-based parameter priors, a model-redundancy prior, and a Bayesian information criterion motivated model prior, penalize complex models to enforce Occam's razor. Using a precomputed grid of simulations, the probabilistic model selection process enables an initial guess for the maximum a posteriori (MAP) material parameter values. Synthetic tests recover the ground-truth models and expected parameters. Using experimental data for gelatin, fibrin, polyacrylamide, and agarose, MAP simulations of credible models reproduce the data. Moreover, a cross-institutional comparison of 10% gelatin indicates consistent constitutive model selection.

---

1 Introduction

Soft biological and tissue-like materials experience extreme mechanical loading during events such as blast exposure, cavitation, and histotripsy, where deformation occurs on the microsecond timescale.^1–4 The combined rapid energy transfer and large deformations during these events place the material response in the high-strain-rate regime (10³–10⁸ s⁻¹).^5–8 Under such conditions, these materials exhibit nonlinear viscoelastic behavior, which is significantly different from their quasistatic response.^9–18 Accurately capturing and characterizing this behavior is important for applications ranging from impact biomechanics to focused ultrasound and soft matter engineering. Yet, it remains difficult due to limited optical access, sensitivity to boundary conditions, and coupled thermal-diffusive effects that arise during rapid loading.^{7,11,12,19–21} Characterization is typically conducted in two parts: the determination of constitutive models and the estimation of parameters. If a constitutive material model has been identified with a priori information in the parameter domain of interest, then this process is simplified to calibrating parameter values.

Conventional experimental approaches only offer partial insight into this extreme regime, but they are limited in several notable ways. Split-Hopkinson pressure bar or direct high-speed impact testing, for example, measures material response at strain rates of (10²–10⁴ s⁻¹)^22–27 but is constrained by inertial effects, weak transmitted signals, and the assumption of stress equilibrium.^28–30 High-speed imaging with digital image correlation (DIC) or particle tracking has been used to capture full-field deformation of soft materials during impact or shock loading.^31–39 However, DIC requires the sample surface to have dense speckle patterns and faces an inherent tradeoff between noise reduction and spatial resolution.^40,41 Particle tracking may introduce mismatch errors.⁴² Computational modeling techniques, such as finite element analysis (FEA) and multiscale constitutive modeling,^43,44 have been employed to simulate the rate-dependent and large-strain behavior of soft materials. Hrapko et al.⁴³ used FEA to model the dynamic behavior of brain tissue using various viscoelastic constitutive models and validated their simulations against shear experiments (0.1–1 s⁻¹). Their analysis revealed a strong sensitivity to model form and parameter values, posing a challenge to achieving consistent and physically reliable model selection in the quasistatic regime. More recently, data-driven and machine learning approaches have been proposed to integrate experimental data into predictive frameworks.^45,46 Although these methods can capture nonlinear and rate-dependent behavior, they often require extensive, high-quality datasets that are rarely available for soft materials.⁴⁷ When data are sparse relative to the model complexity, overfitting can occur, while limited coverage of the relevant parameter or strain-rate space can lead to underfitting.⁴⁸

Cavitation has proven to be an effective rheometry technique for probing the mechanical response of soft materials across a wide range of strain rates.¹³ Early work by Zimberlin et al.⁴⁹ showed that needle-induced cavitation in hydrogels can be used to infer local stiffness. Subsequent studies extended these approaches to biological tissues, including brain tissue, where cavitation dynamics have been used to extract viscoelastic and failure properties.^5,50 More recent work has shown that inertial cavitation can provide meaningful rheological information in hydrogels at high strain rates. For example, Cohen and Durban⁵¹ analyzed the dynamic expansion of a spherical cavity in compressible elastoplastic materials and identified conditions under which plastic shock waves appear. Related efforts in soft viscoelastic media by Warnez and Johnsen,⁵² Gaudron et al.⁵³ incorporated rate-dependent dissipation and elastic nonlinearity into cavitation dynamics models, showing that viscosity can delay rapid cavity growth and attenuate shock-like responses.

Leveraging these prior insights, the Inertial Microcavitation Rheometry (IMR) technique was developed to enable the characterization of soft hydrogels across strain rates ranging from 10³–10⁸ s⁻¹.⁶ Under varying pressures and material constraints,^{8,19,20,54–58} the radial dynamics of microcavitation are modeled using modified Keller–Miksis equations. In IMR, materials are characterized by comparing laser-induced microcavitation experimental observations to Keller–Miksis-type forward numerical calculations for simple viscoelastic constitutive models using a least-squares fitting approach.⁶ Experimental radius–time data are typically extracted near the rebound peaks, when the bubble wall velocity is low and optical measurements are most reliable. High-speed imaging at frame rates greater than one million frames per second captures most of the bubble dynamics. However, rapid motion near collapse may fall outside the temporal resolution or show motion blur.

The IMR technique considers multiple viscoelastic constitutive material models^7,19,57,59 to characterize hydrogels, e.g., gelatin, agarose, fibrin, and polyacrylamide. Spratt et al.⁶⁰ extended the technique using ensemble-based data assimilation methods, obtaining similar mechanical properties for the same materials as Estrada et al.⁶ while reducing computational cost. Their results suggested that achieving a best fit between simulation and experiment may require a viscosity that evolves during the bubble collapse, potentially reflecting damage or structural softening effects not captured by conventional constitutive models.⁶⁰ Zhu et al.⁵⁸ developed a parsimonious IMR technique that simplifies soft material characterization by focusing on the initial bubble collapse time rather than fitting the entire radius–time history. While these advances improved the IMR methodology, consistently identifying the best constitutive model capturing the material response at these high strain rates remains a challenge. In practice, candidate models are evaluated sequentially using heuristic decision procedures with least-squares fitting. Although these approaches can produce reasonable parameter estimates for individual models, they do not quantitatively compare the models or assess uncertainty associated with the model selection.

Bayesian inference provides a probabilistic approach for extending the IMR technique to model selection and parameter estimation. Unlike deterministic least-squares methods that yield a single best-fit solution, Bayesian inference can quantify the plausibility of competing constitutive models by evaluating model performances through their relative probabilities^61–64 and then identify the associated parameters.^65–69 This probabilistic formulation directly evaluates model accuracy and quantifies uncertainty without relying on extensive prior information.^70–72 Additionally, Bayesian inference naturally incorporates model-form priors, experimental uncertainty, and measurement noise across multiple datasets,^71,73,74 making it well-suited for reproducible and uncertainty-aware characterization of soft materials.

Recent studies have shown the growing impact of Bayesian inference in rheological problems. Chu et al.^75,76 applied Bayesian optimal experimental design to bubble dynamics and soft material characterization, showing that information-theoretic criteria can accelerate data acquisition and improve inference efficiency. Similarly, Freund and Ewoldt⁷⁷ used Bayesian model selection to determine the optimal number of modes in a multi-modal Maxwell model, which describes the dynamic shear moduli of a synthetic polymer. Their results showed that Bayesian inference can reliably discriminate between competing viscoelastic models with high-dimensional parameter spaces and noisy experimental data. In this context, Bayesian inference favors constitutive models that strike a balance between goodness of fit and physical interpretability, rather than empirical or highly flexible formulations with many parameters that risk overfitting. While effective in selecting a model, their method was unable to quantify uncertainty. Further work extended Bayesian approaches to rheological characterization, including hierarchical Bayesian inference for viscoelastic parameter estimation and Bayesian uncertainty quantification for constitutive modeling.^78,79 Data-driven frameworks have also explored the identification of rheological constitutive behavior from partial flow measurements and differentiable simulations, enabling model discovery from velocity or stress observations in complex flows.⁸⁰ However, these approaches primarily consider experimental configurations in which stress–strain or flow-field measurements are observed directly. On the other hand, IMR experiments only provide transient bubble radius–time dynamics, requiring material properties and constitutive behavior to be inferred indirectly from the bubble dynamics rather than from directly measured stress or strain fields.

To this end, this work aims to implement Bayesian inference in conjunction with the IMR technique to perform model selection across competing constitutive models under high-strain-rate conditions, while accounting for uncertainty quantification. The work is organized as follows. In Section 2, we present the governing equations for the bubble dynamics used to generate forward simulations across a parameter space of constitutive model parameters, the experimental setup for laser-induced cavitation (LIC), the material characterization procedure, and the Bayesian inference IMR implementation. The results of the Bayesian IMR characterization for soft materials and discussion of their implications are presented in Sections 3 and 4, respectively. We summarize the capabilities of the Bayesian IMR technique and outlook in Section 5.

2 Theory and methods

2.1 Governing equations

We consider a spherical bubble that is nucleated in a soft material by a laser-induced cavitation event. As shown in Fig. 1, the bubble expands to a maximum radius and subsequently collapses and oscillates toward a mechanical equilibrium. To model these dynamics, the IMR forward simulation code base is used to perform numerical simulations of the governing equations that describe bubble oscillations. The equations consist of an ordinary differential equation (ODE) for the bubble's radial dynamics, which is coupled with partial differential equations (PDEs) for the energy balance inside the bubble, and a PDE for mass transfer between the bubble's water vapor and non-condensible gas.^21,52 A numerical resolution of 100 points inside the bubble was used for the forward simulations. In the present study, stress evolves through explicit functions or ODEs, as defined in Table 1.^6,53,58 For brevity, only the radial equation for the bubble is stated here. The bubble's radial dynamics are modeled using the Keller–Miksis equation:^81,82where the dot denotes a time derivative, R is the bubble radius, Ṙ the bubble wall velocity, c the speed of sound in the surrounding material, ρ its density, S the material viscoelastic contribution of the stress integral, γ the surface tension between the material and the vapor and gas filled bubble, and p_b and p_∞ the bubble and far-field pressures, respectively. S is defined aswhere τ_rr and τ_θθ are the respective radial and hoop components of the deviatoric Cauchy stress predicted by the chosen constitutive model, and r the radial coordinate in the surrounding medium.

Fig. 1 A spherical bubble is nucleated in a soft material and grows to a maximum radius before collapsing and oscillating until it reaches mechanical equilibrium. Identification of the stress integral formulation enables the characterization of the material.

Table 1 Constitutive material models and corresponding stress integral

Tag	Constitutive model
Newt	Newtonian fluid
NH	Neo-Hookean solid
KV	Kelvin–Voigt solid
qNH	Quadratic neo-Hookean solid
LM	Linear Maxwell fluid
qKV	Quadratic Kelvin–Voigt solid
SLS	Standard linear solid

---

For comparison across materials, material models, and experimental conditions, it is beneficial to non-dimensionalize the governing equations. We present the governing equations in non-dimensional form. We choose the bubble radius at its maximum size, L_c = R_max, as the characteristic length, the far-field pressure and material density to define the characteristic velocity, , and the characteristic time, .^6,55,83 The non-dimensionalized variables are

Thus, the non-dimensional Keller–Miksis equation,⁶

The explicit forms of S* for each constitutive model are tabulated in Table 1, expressed in terms of the wall stretch ratio , where denotes the non-dimensional equilibrium bubble radius. The dimensionless parameters are the Weber number, Reynolds number, Cauchy number, and Deborah number defined as,

respectively, where μ is the surrounding material shear viscosity, G the elastic shear modulus, λ₁ = μ/G₁ is the relaxation time. The strain-stiffening parameter, α, is also non-dimensional. The constitutive model parameter ranges are tabulated in Table 2. Within these bounds, the parameter space is sampled in log space to ensure uniform coverage across orders of magnitude for the simulation campaign. We fix the total number of forward simulations per model to 4096 and construct a Cartesian grid with n^d = 4096 points, where d is the number of model parameters and n = 4096^1/d. Consequently, the per-axis resolution decreases with increasing dimensionality (e.g., n = 4096 for d = 1, n = 64 for d = 2, n = 16 for d = 3), ensuring comparable total computational effort across models of different dimensions. A convergence study of the observational likelihood with respect to grid resolution is provided in Section A of the appendix, demonstrating that the computed likelihood values stabilize at the chosen resolution.

Table 2 Simulation material parameters for soft-tissues in the high-strain rates and their non-dimensional counterparts^6,15,60

Parameter	Notation	Minimum	Maximum
Viscosity [Pa s]	μ	10^−4	1
Reynolds number	Re	1.6	36 × 10^3
Shear modulus [Pa]	G	10^2	10^5
Cauchy number	Ca	0.2	10^3
Relaxation time [s]	λ1	10^−7	10^−3
Deborah number	De	3 × 10^−3	65
Strain stiffening	α	10^−3	10

---

2.2 Experiments

LIC experiments in soft materials were conducted and captured using a high-speed camera.^4,7,8,84 The University of Michigan (UM) and the University of Texas at Austin (UT) LIC experimental setups, material preparation protocols, and datasets for the hydrogels are in ref. 57, 59 and 85, respectively. The materials and their corresponding polymer concentrations are tabulated in Table 3, including chemically crosslinked polyacrylamide (PAAm; 8 wt% acrylamide with 0.26 wt% bisacrylamide) and physically crosslinked gelatin (10 wt%), fibrin (0.2 wt%), and agarose (5 wt%). The cavitation experiments were performed as single, spatially isolated events to avoid damage accumulation. Trials exhibiting visible macroscopic fracture, loss of optical tracking, or non-recovering behavior were excluded during experimental preprocessing. A comparison of the maximum bubble radius R_max with the maximum stretch ratio Λ = max(λ) is shown in Fig. 2. Each dataset consists of multiple experimental trials for one material. A subset of UM hydrogel concentrations were used in the present work. Gelatin data of the same concentration from both institutions were analyzed to cross-validate the constitutive model selection. The experimental setup generates and records a single LIC event in transparent hydrogels at ultra-high-speed frame rates (≥10⁶ fps). The image processing is done using an in-house bubble edge detection MATLAB code.⁸⁶

Table 3 Institution, materials, and concentrations of the laser-induced microcavitation experiments

Tag	Material	Concentration [wt%]
UM1	Gelatin	10.0
UM2	Fibrin	0.2
UM3	PAAm	8.0/0.26
UT1	Gelatin	10.0
UT2	Agarose	5.0

---

Fig. 2 Experimental maximum stretch ratio vs. maximum radius. UM1: blue circle, UM2: red square, UM3: black up triangle, UT1: cyan right triangle, UT2: green down triangle.

In Fig. 3, the left panel shows normalized bubble radius histories for PAAm, where differences between trials primarily reflect variations in the applied laser energy, with added contributions from measurement noise, slight material heterogeneity, and minor differences in viscous response due to varying R_max. The middle panel shows the corresponding bubble wall velocities, computed from the non-dimensional experimental R* data using a Padé derivative approximation. The bubble wall radius and velocities are used to compute the non-dimensional strain rate and strain via the logarithmic Hencky strain defined in eqn (4). The right panel shows the phase space between the non-dimensional strain rate and strain, which are considered for filtering out low-information data points in Section 2.2.1. Although the normalized LIC curves collapse into a consistent temporal window, variations in the collapse and rebound amplitudes across trials show the combined effects of experimental noise and material variability. These observations motivate the use of probabilistic model comparison, as deterministic fits do not adequately represent the observed data spread.

Fig. 3 Bubble radius (left), velocity (middle), and magnitude of strain rates at the bubble wall (right) histories for nine sets of PAAm LIC experimental data.

2.2.1 Pre-processing. For a given soft material, there is a dataset of J experimental trials (radius–time curves) indexed as j = 1,…,J. Each trial is non-dimensionalized by its own peak radius, R_max,j, and characteristic time, t_c,j. After pre-processing, each trials provides a set of non-dimensional time stamps and corresponding bubble radius and bubble wall velocity measurements , where and [scr K, script letter K]_j denotes the index set of valid time-step indices for trial j. In some trials, the experimental signal extends far beyond the collapse events of interest, with small oscillations about equilibrium that contribute disproportionately to the Gaussian likelihood. To focus the Bayesian analysis on informative regions, [scr K, script letter K]_j is restricted using a combination of strain-based filters. We define the radial normal component of both the Hencky strain with respect to the per-trial equilibrium radius, , and the dimensionless strain rate as

We set threshold values at

where the strain threshold corresponds to the onset of finite-strain behavior in soft hydrogels, beyond which small-strain assumptions no longer hold,⁸⁷ and the strain rate threshold corresponds to the high-strain-rate regime. Time-steps are retained if they lie outside an elliptical gate in the plane,which removes near-equilibrium points characterized by low strain and low strain rate, while preserving the collapse-rebound regions of interest.

Each simulation is generated using the mean values extracted from the experiments,

For a given model M_i and associated parameter vector θ_i, where i ∈ [1,N_M] and N_M is the number of models, the forward simulation produces a discrete time history , where the subscript s indicates simulation data and m = 1,…,N_t indexes the simulation time steps with maximum N_t time steps. The simulation data are then interpolated to the experimental time stamps to obtain the simulated radius and velocity at the proper time steps, denoted and . Residuals between the experimental observations and the simulated predictions are defined as

where θ_i = {θ₁,…,θ_n} is the set of n model parameters for model M_i.

2.3 Characterization

For a single material and fixed laser input energy, repeated LIC events reveal variability between experimental realizations. A 95% confidence interval computed across multiple trials provides statistical bounds on the expected variation at each time step. However, the underlying noiseless response may not be fully captured within the available experimental datasets. While the pre-processing step filters out low-information regions to limit their influence, residual uncertainty remains in both the retained and discarded portions of the data. The Bayesian inference approach introduced in Section 2.3.1 incorporates this uncertainty through the likelihood formulation, which explicitly accounts for observational noise and potential model–data mismatch. Since the forward simulations solve deterministic governing equations, the modeling noise is assumed to be zero. However, the probabilistic formulation still enables consistent treatment of measurement variability across trials.

2.3.1 Bayes' theorem for model selection. A Bayesian approach is used to select a constitutive model and quantify the agreement between the IMR forward simulation results and the LIC experimental data. Bayes' theorem for model selection is^88–90where P(M_i|D,I) is the posterior probability of a model, P(M_i|I) the model prior, P(D|M_i,I) the observational likelihood, and P(D|I) the marginal likelihood for i-th model, M_i, D the given data, and I the system constraints. For simplicity and consistency with the literature, the constraint term is omitted from the notation, though it is understood to be present. Additionally, an intrinsic Occam's razor penalizes overly complex models, yielding the simplest model that maximizes the observational likelihood, where complexity refers to the number of model parameters.⁹¹ For IMR characterization, D is the experimental and numerical simulation bubble radius–time and bubble wall velocity-time data, i.e., R*(t*) and Ṙ*(t*).

2.3.2 Hierarchical noise scaling and likelihood. Baseline heteroscedastic variance estimates are computed per time step, denoted by and . To model the time-varying reliability of experimental measurements, we apply a logistic weight that depends on the instantaneous strain rate,where κ controls the steepness of the transition around the threshold strain rate . κ was set to unity to have a gradual transition that is neither abrupt nor overly diffuse, allowing the weighting to change on a physically meaningful scale comparable to typical variations in measured strain rate. Tests with larger κ values produced sharper transitions that offered no improvement in model discrimination and occasionally reduced numerical robustness. The logistic function is then mapped to a bounded weight w_j,k via

w_j,k = m + (1 − m)a_j,k, m ∈ (0,1], κ > 0, (11)

where m sets the minimum allowable weight, with the maximum being unity. The value of m is set to 0.1, limiting down-weighting for data points in low-activity regions by at most a factor of ten relative to the highest-weighted points. This small floor prevents near-zero weights from inflating the variance. In this formulation, measurement uncertainty in high-speed optical imaging increases primarily with bubble wall velocity, and thus with strain rate. Slow rebound or near-equilibrium dynamics yield minimal blur and tracking error, while rapid bubble collapses result in reduced resolution due to motion and finite camera exposure. Therefore, including strain amplitude in the weighting would penalize optically reliable, slower-moving segments and misrepresent the true measurement noise. In contrast, strain-rate-based weighting isolates the physically dominant source of heteroscedasticity without redundantly penalizing regions already excluded by the elliptical gate defined above. The scaled variances arewhere β is a positive noise scale such that β > 0 and inflates the per-time experimental variance. The inferred noise scale β is not a direct estimate of experimental measurement noise. It is an effective discrepancy parameter absorbing mismatch sources between the model and the data (e.g., measurement noise, under-resolved dynamics, and model inadequacy.) Furthermore, β serves as a diagnostic of data quality and as a safeguard against overfitting. Models that require artificially large β values to reconcile simulations with data are automatically downweighted in the plausibility calculation. Thus, as shown in eqn (15), high posterior probabilities are reserved for models that explain the data within a realistic noise margin. This interpretability is particularly valuable when multiple viscoelastic formulations achieve similar log-likelihoods. The combined consideration of plausibility and β distributions distinguishes between models that genuinely represent the dynamics and those that rely on the inflated noise scaling.

A simulation/model-discrepancy variance, and , can be included to account for uncertainty in stochastic or approximate forward models. The total variances are then

Here, the simulation/model-discrepancy variance terms are set to zero by default because the IMR forward simulations are deterministic.

To obtain P(D|M_i), the parameters and noise scale are marginalized. When modeling physical systems, Gaussian noise is commonly introduced to account for the inherent experimental variability and uncertainty, maximizing statistical entropy.^88,92 A Gaussian observational likelihood is adopted,

2.3.3 Priors. We employ a hierarchical prior structure that acts at three coupled levels: the noise-scale prior P(β), the parameter prior P(θ_i|M_i) within each constitutive model, and the model prior P(M_i) across the set of candidate models. In principle, one could adopt a uniform prior over a considered space, P(·) = 1/V, where V is either the total number of elements or the volume in the space, depending on if it is discrete or continuous. Such a prior corresponds to the state of maximum statistical entropy and is therefore the least biased choice in the absence of prior information.^70,88,92,93 However, when models differ in dimensionality or physical formulation, the uniform assumption may not sufficiently penalize model complexity and assign similar plausibility to parameter regions that are physically redundant or poorly constrained. Therefore, we construct priors to penalize model complexity. The noise prior constrains the level of variance inflation permitted when reconciling simulations with experimental data. The parameter space prior ensures consistent weighting across physically meaningful and unique regions of the parameter space. The model prior then imposes an Occam penalty to prevent overparameterization.

To marginalize out β from eqn (14), we compute

whereis the half-Cauchy distribution. A half-Cauchy prior for β is a weakly informative, heavy-tailed distribution that regularizes small scales while allowing large values when supported by the data. Thus, it is a robust default for variance components in hierarchical Bayesian models.⁹⁴

Because the integral in eqn (15) has no closed form antiderivative for the likelihood function used here, it is evaluated numerically by discretizing the domain of β into B grid points {β_b}^B_b=1,

where the normalized quadrature weights are defined asand Δβ is the uniform grid spacing. Although the half-Cauchy prior is defined on a semi-infinite domain, its heavy-tailed form decays as P(β) ∼ β⁻² for large β, meaning that most of its probability mass is concentrated near β ≲ 10. Thus, we truncate the integration domain to the finite interval [β_min, β_max] = [0.05, 10], which captures more than 99.9% of the total prior mass with negligible truncation error. Within this finite range, the quadrature is performed on a uniform grid using the trapezoidal rule, with each grid point being the midpoint of a local interval between its neighbors. This symmetric construction preserves second-order accuracy and ensures that both interior and boundary points contribute proportionally to the prior-weighted integral.

To marginalize out θ_i, we compute

Similar to eqn (15), P(D|M_i,θ_i) has no closed-form antiderivative in θ_i for the likelihoods considered here. Thus, we evaluate eqn (19) numerically on a Cartesian grid :

where g indexes each unique grid point and N_g is the total number of parameter combinations. The quadrature weights w_g represent the normalized prior mass assigned to the cell at θ^(g)_i:where ΔV_g is the (hyper-)volume of the grid cell.

We now define the prior for the parameter space to compute eqn (21). For each constitutive model M_i, the parameter grid θ_i is mapped to a normalized space using logarithmic coordinates. This normalization makes the parameter representation range-invariant with respect to the choice of parameter bounds. The model's stress integral response is evaluated along the gated experimental segments defined in eqn (6), denoted by . These stress integrals are used to determine whether additional parameters of M_i generate distinct stress integral histories relative to lower-dimensional models in the hierarchy shown in Fig. 4.

Fig. 4 Schematic diagram of the hierarchical relationship among the constitutive models. Arrows indicate limiting cases obtained by setting parameters to zero, thereby reducing the model's dimensionality. The first column of models corresponds to three parameters, the second column has two parameters, and the third column has one parameter.

To quantify this distinction, we consider parent and child models, M_{[scr P, script letter P]} and M_{[capital script C]}, respectively, such that M_{[capital script C]} ⊆ M_{[scr P, script letter P]} and θ_{[capital script C]} ⊆ θ_{[scr P, script letter P]} represents the corresponding parameter subsets. Models are distinguished by differences in their respective stress values. Therefore, we consider the relative stress mismatch between two models,

where ‖·‖_w denotes the weighted ℓ² norm induced by {w_j,k}. Small values of Φ_{[scr P, script letter P]→[capital script C]} indicate that the additional parameters in M_{[scr P, script letter P]} do not introduce stress integral histories that are distinguishable from those of the simpler model M_{[capital script C]}.

To convert this mismatch into a prior penalty, we define a dimensionless redundancy factor

where τ_{[capital script C]} is a data-driven stress scale defined from the median and median–absolute–deviation of the simulated stress signal, representing the smallest stress change the model can reliably resolve. Differences smaller than τ_{[capital script C]} are not resolvable by M_{[capital script C]} and therefore do not constitute genuine distinctions between the two models. To ensure that M_{[scr P, script letter P]} is penalized whenever it can emulate any lower-dimensional model, we take

Finally, the prior for the parameter space is constructed by combining a harmonic-mean bottleneck of the form

where d is the model dimension and a small, machine precision value of ε > 0 is included for numerical stability, and the redundancy penalty,

P(θ_i|M_i) = H_i(θ_i)w_red(θ_i), (26)

such that . The sensitivity of the redundancy-based parameter prior to reasonable variations in the redundancy definition is examined in Section B of the appendix.

Finally, we define the model prior. Similar to the parameter prior, we consider the dimensionality of the models and their effect on the posterior. We construct a weakly informative prior that penalizes model complexity,

where k_M is the number of free parameters in model M_i and N_eff is the effective number of scalar observations. This prior is analogous to the Bayesian information criterion complexity term and introduces an explicit Occam penalty.⁹⁵

2.3.4 Posterior. To compute the posterior P(M_i|D), we consider the marginal likelihood. However, without additional prior information about the data, calculating the marginal likelihood is a challenging task. Therefore, for the model comparison process, P(D) is constant in the models to normalize the posteriors and, thus, the sum of the posteriors is unity:⁸⁹

The likelihood and prior values can span orders of magnitude across grid points and trials, making direct summation numerically unstable. To maintain numerical precision and prevent underflow, multiplicative terms are computed in log space, and marginalizations are carried out using the log-sum-exp operation. Taking the logarithm of eqn (9) yields

where log P(M_i) penalizes models according to prior information, log P(D|M_i) rewards models that fit the data well, and log P(D) serves as the normalization constant ensuring that posterior probabilities sum to unity. The effect of Occam's razor becomes explicit in log space: the prior term penalizes models with greater complexity (i.e., larger parameter dimension k_M = d) relative to the effective number of observations N_eff.

2.3.5 Additional quantities of interest. After selecting the most plausible model M*, we compute two additional posterior distributions of interest using the same likelihood evaluations and parameter grid employed during the computation of the model likelihood P(D|M_i). The observational likelihood, quantifies the probability of observing the data under a given model after integrating out both the model parameters and the noise scale:

This integral is approximated using the same discrete grids in θ_i and β that were used for the likelihood evaluations, so no additional forward simulations are required.

The marginal posterior over parameters θ_g is obtained by integrating out the noise scale β using normalized prior weights over the discrete β grid {β_b}^B_b=1:^70,88

where the normalized prior weights w_b approximate the prior mass over β via

These weights arise from discretizing the continuous marginalization integral over β (eqn (15)) into a quadrature sum, where each w_b corresponds to the probability content of the local interval around β_b under the prior P(β).

The prior over θ_g is uniform across grid points. The resulting discrete posterior distribution is normalized such that

and the maximum a posteriori (MAP) parameter estimate, θ_MAP, corresponds to the grid point g maximizing eqn (28).

The posterior over β for a fixed parameter value θ_g is given by

P(β_b|D,M_i,θ_g) ∝ P(D|M_i,θ_g,β_b)P(β_b), (29)

which is normalized over b so that . When θ_g is chosen as θ_MAP, the MAP estimate β_MAP is simply the β_b that maximizes eqn (29). Maximizing the posterior over β provides a principled estimate of the most likely noise scale in the data for the chosen model, which serves as a direct measure for the uncertainty quantification (UQ) of the fitted dynamics. A larger β_MAP represents greater experimental variability and correspondingly wider predictive intervals. In contrast, a smaller β_MAP reflects higher confidence in the model's ability to reproduce the observed data within the inferred noise bounds.

Algorithm 1 shows the hierarchical Bayesian inference workflow, including heteroscedastic likelihood evaluation, marginalization over the noise scale parameter β, model likelihood computation, and extraction of θ_MAP and β_MAP. The algorithm also outputs the posterior distributions P(β|D,M_i) and P(θ_i|D,M_i) for UQ.

Algorithm 1 Hierarchical Bayesian model selection for IMR

Input: Experimental data ; set of candidate models M_i with parameter grid {θ_i,g}; uniform noise–scale grid {β_b} with prior P(β)

(1) Compute baseline variances and ; strain rates ; activity a_j,k; and weights w_j,k

(2) Using the stress–integral mismatch measures Φ_{[scr P, script letter P]→[capital script C]} and redundancy factors w_{[scr P, script letter P]→[capital script C]} (see eqn (22) and (23)), compute the parameter–space priors P(θ_i|M_i) for all models via eqn (26)

(3) for each model M_i do

(4) for each grid point θ_i,g do

(5)  Interpolate simulations to trial j's time stamps

(6)  for each β_b do

(7)   Evaluate the observational likelihood P(D|M_i,θ_i,g,β_b) using the Gaussian likelihood with heteroscedastic variance (see eqn (14))

(8)  end for

(9)  Marginalize over β using eqn (15) to obtain P(D|M_i,θ_i,g)

(10) end for

(11) Compute the model likelihood P(D|M_i) via eqn (20) by summing over θ_i,g

(12) end for

(13) Compute posterior model probabilities P(M_i|D) from Bayes' theorem eqn (9) and select the most plausible model M*

(14) Optionally compute the marginal parameter posterior P(θ_i,g|D,M_i) (see eqn (28)) and the noise–scale posterior P(β|D,M_i) (see eqn (29)) to obtain θ_MAP and β_MAP

Output: P(M_i|D), P(θ_i,g|D,M_i), P(β|D,M_i), selected model M*, θ_MAP, β_MAP

3 Results

3.1 Consistency check

We perform a consistency check to recover the model and parameter values of synthetic data for which the true viscoelastic model and parameters are known a priori. The consistent model selection from synthetic ground-truth conditions ensure high fidelity inverse characterization of the experimental data. Using the bubble dynamics model described in Section 2.1 along with the stress integral formulations in Table 1, 32 sets of R*(t*) synthetic data for each of the constitutive models are simulated. The parameter values used to generate synthetic data were chosen for their relevance to the soft materials of interest. Specifically, the values used were μ = 0.05 Pa s, G = 10 kPa, λ₁ = 10⁻⁵ s, and α = 1. R_max is varied with respect to the observed experimental LIC data from Fig. 2. Each constitutive model is then compared to the synthetic data using the approach in Section 2.3. The results are tabulated in Table 4. The method correctly identifies generative models that serve as the ground-truth with close parameter values and posterior values equal to unity.

Table 4 Consistency check using synthetic data with known parameter values. The rows highlighted in bold represent the most plausible model with a posterior value of unity

Ground-truth	Model	μ (Pa s)	G (kPa)	λ1 (s)	α
	Newt	0.051	—	—	—
	NH	—	10.6	—	—
Newt	KV	0.046	0.58	—	—
μ = 0.05 Pa s	qNH	—	9.91	—	0.009
	LM	0.072	—	1.0 × 10^−7	—
	qKV	0.046	0.55	—	0.012
	SLS	0.046	3.02	1.0 × 10^−7	—
	Newt	0.049	—	—	—
	NH	—	10.4	—	—
NH	KV	0.001	9.91	—	—
G = 10 kPa	qNH	—	9.91	—	0.014
	LM	0.062	—	1.0 × 10^−7	—
	qKV	0.007	3.02	—	0.464
	SLS	0.014	9.39	1.85 × 10^−7	—
	Newt	0.138	—	—	—
	NH	—	20.4	—	—
KV	KV	0.054	9.91	—	—
μ = 0.05 Pa s	qNH	—	19.5	—	0.001
G = 10 kPa	LM	0.173	—	1.79 × 10^−7	—
	qKV	0.046	9.39	—	0.022
	SLS	0.086	9.39	1.0 × 10^−7	—
	Newt	0.423	—	—	—
	NH	—	88.0	—	—
qNH	KV	0.0193	65.8	—	—
G = 10 kPa	qNH	—	13.0	—	0.720
α = 1	LM	0.268	—	8.96 × 10^−7	—
	qKV	0.001	16.6	—	0.464
	SLS	0.293	91.0	1.17 × 10^−6	—
	Newt	0.001	—	—	—
	NH	—	5.54	—	—
LM	KV	0.001	5.77	—	—
μ = 0.05 Pa s	qNH	—	4.41	—	0.093
λ1 = 10^−5 s	LM	0.054	—	5.99 × 10^−6	—
	qKV	0.001	3.02	—	0.251
	SLS	0.086	1.71	1.36 × 10^−5	—
	Newt	0.423	—	—	—
	NH	—	57.1	—	—
qKV	KV	0.072	75.3	—	—
μ = 0.05 Pa s	qNH	—	38.3	—	0.125
G = 10 kPa	LM	0.359	—	2.78 × 10^−7	—
α = 1	qKV	0.046	16.6	—	0.464
	SLS	0.046	51.6	1.0 × 10^−7	—
	Newt	0.373	—	—	—
	NH	—	2.56	—	—
SLS	KV	0.481	29.2	—	—
μ = 0.05 Pa s	qNH	—	7.57	—	1.49
G = 10 kPa	LM	0.359	—	1.0 × 10^−7	—
λ1 = 10^−5 s	qKV	0.541	29.2	—	0.002
	SLS	0.293	29.2	1.36 × 10^−5	—

---

3.2 Inverse characterization via Bayesian model selection

The Bayesian inference method was applied to the experimental UM and UT LIC datasets. The model plausibilities and MAP parameter estimates are tabulated in Table 5. The optimal noise scaling parameter is included for the experimental data results to account for the variability of the experimental data.

Table 5 Plausibilities for UM and UT experimental data and corresponding parameter values for each model. The rows highlighted in bold represent the most plausible model

Material	Model	μ (Pa s)	G (kPa)	λ1 (s)	α
UM1	Newt	0.058	—	—	—
	NH	—	20.2	—	—
	KV	0.022	13.0	—	—
	qNH	—	19.5	—	0.01
	LM	0.083	—	1.16 × 10^−7	—
	qKV	0.025	9.39	—	0.074
	SLS	0.046	16.6	1.85 × 10^−7	—
UM2	Newt	0.148	—	—	—
	NH	—	7.71	—	—
	KV	0.173	9.91	—	—
	qNH	—	7.57	—	0.003
	LM	0.173	—	1.0 × 10^−7	—
	qKV	0.158	9.39	—	0.003
	SLS	0.001	1.71	4.64 × 10^−5	—
UM3	Newt	0.195	—	—	—
	NH	—	41.6	—	—
	KV	0.112	22.3	—	—
	qNH	—	38.3	—	0.007
	LM	0.173	—	1.16 × 10^−7	—
	qKV	0.158	5.32	—	0.464
	SLS	0.086	29.2	1.0 × 10^−7	—
UT1	Newt	0.125	—	—	—
	NH	—	45.2	—	—
	KV	0.072	33.5	—	—
	qNH	—	43.9	—	0.012
	LM	0.645	—	6.94 × 10^−6	—
	qKV	0.086	29.2	—	0.074
	SLS	0.086	29.2	1.0 × 10^−7	—
UT2	Newt	0.441	—	—	—
	NH	—	166.7	—	—
	KV	0.645	381.5	—	—
	qNH	—	13.0	—	4.81
	LM	0.416	—	2.40 × 10^−7	—
	qKV	0.158	9.39	—	10.0
	SLS	1.00	500	3.41 × 10^−7	—

---

Most datasets consistently placed the posterior mass on the KV model, with a single dataset favoring the qKV model. For UM1, KV had a posterior of 0.83 and qKV had a posterior of 0.17. For the other cases, the highlighted model in Table 5 had a posterior of unity. The finding that multiple materials select the same constitutive family does not imply identical behavior. Their MAP parameter estimates occupy distinct and well-separated regions of parameter space, reflecting genuine material differences. In particular, the KV-selected materials differed primarily in their relative elastic-to-viscous values, shifting collapse-time sensitivity and rebound curvature in characteristic ways, as shown in the left column of Fig. 5. The dataset selecting qKV showed a sharper rebound curvature and steeper post-collapse trajectory compared to the others, which the linear KV model could not reproduce. The Bayesian likelihood correctly identified this mismatch and shifted posterior mass toward the nonlinear model.

Fig. 5 Bayesian IMR characterization for experimental datasets. Left column: Normalized bubble radius histories R* versus non-dimensional time t* showing the experimental data cloud (cyan), the MAP simulation of the most plausible model (red), and the posterior predictive band (purple) representing the 5–95 percentile interval obtained from nearby parameter realizations. Middle column: Normalized absolute residuals showing the time-dependent discrepancy between experimental observations and the MAP prediction. Right column: Posterior uncertainty quantification for the noise scale parameter β in the selected model M*, with the marginal P(β|D,M*) (solid blue) and the conditional posterior θ_MAP (dashed red).

To evaluate the quality of the MAP estimates, we compute posterior predictive intervals and assess their empirical coverage. The latter is the fraction of experimental observations contained within the 5–95 percentile posterior predictive band. The band is constructed by collecting precomputed simulations at ±2 discrete points of θ_MAP and injecting heteroscedastic noise at each time step. Since the parameter grid is sampled in log space, each step corresponds to a fixed increment, resulting in the neighborhood spanning a physically meaningful range around the MAP. Coverage values are 96.1%, 94.7%, and 90.0% for UM1–UM3 and 100.0% and 87.9% for UT1–UT2, indicating that the inferred models reproduce the observed trial variability across datasets. The middle column of Fig. 5 shows that the absolute residuals remain small for most of the bubble evolution, with larger deviations appearing during collapse and rebound phases.

Finally, we examined the inferred noise-scale β posteriors across materials. Rather than acting as a free parameter that compensates for model mismatch, β reflects how consistently the experimental trajectories align with the model-predicted collapse-rebound dynamics. Datasets with trials following similar radius histories produce a sharply concentrated posterior for β, indicating that only a modest amount of measurement noise is required to explain the observed variation. In contrast, datasets with visibly greater scatter yield broader noise-scale posteriors, reflecting genuine variability in the recorded bubble motion. The right column of Fig. 5 compares the marginal noise-scale posterior (solid blue), P(β|D,M*), with the conditional posterior at the MAP parameters (dashed red), P(β|D,M*,θ_MAP). With the exception of the UT2 dataset, the red curves are broader and shifted toward larger β than the blue curves, which are narrower and lower in value.

4 Discussion

4.1 Synthetic data insights

The synthetic datasets provide a controlled environment for evaluating the Bayesian IMR method under known ground truths, showing both the strengths and limitations of the technique. Across the synthetic cases, the method correctly recovered the true generating model and appropriate parameter estimates, but only when the informed priors were used. Under uniform priors, the model likelihood is dominated by the volume of the parameter space rather than the quality of fit. In several synthetic tests, this led to incorrect model selection and favored higher-dimensional models (data not shown). The constructed priors resolve this issue by weighting the parameter space according to the relative information. Regions that mimic simpler models receive reduced prior mass, preventing inflated support for unnecessarily complex models. Specifically, the Bayesian likelihood penalizes regions of the parameter space unless they produce a genuine improvement in the dynamics. The synthetic tests thus confirm that the priors are not an optional regularizer but an essential component that aligns the inference with the underlying physics.

Another consideration is the effect of the parameter-grid resolution on the computed likelihoods. Since the likelihood is integrated across the entire grid, models with more parameters require exponentially more grid points to cover their domain with comparable resolution. If the models use the same course grid (e.g., 16 values per dimension), the higher-dimensional models incorrectly dominate (data not shown). This occurs because course grids over-smooth the likelihood landscape of the simpler models while still providing sufficient coverage for the complex ones, artificially inflating their likelihood. For this reason, we use model-dependent, computationally tractable grids.

4.2 Experimental findings

The characterization of the experimental datasets reveals the strengths of the Bayesian IMR technique. The selected model and parameters matched the early-time collapse and rebound envelopes of the experimental clouds, with most materials placing their posterior mass on the KV model. For the cross-institutional gelatin study, the Bayesian IMR technique identified the KV model as well for similar experimental setups and data acquisition methods. The widespread selection of KV is due to its simplicity and accuracy, as additional parameters from the other considered models did not provide any improvement in likelihood. The only exception is the UT2 dataset, where the posterior favors the qKV model. This shift arises from a reproducible feature in the first rebound where the data shows a much sharper recoil than the other materials. The KV model cannot generate this feature without distorting the bubble radius history after the second collapse. The nonlinear strain-stiffening term correctly adjusts the rebound slope without compromising the collapse timing. Overall, these results confirm that the technique extracts constitutive behavior that is physically interpretable under realistic variability.

The posterior predictive intervals further support this interpretation. Their empirical coverage show the inferred models reproducing the observed variability across experimental trials. The remaining discrepancies are concentrated primarily near collapse events, where the bubble dynamics evolve most rapidly. Otherwise, the MAP simulations match experimental trajectories, indicating that inferred constitutive models capture the dominant material viscoelastic response.

4.3 UQ interpretability

The posteriors for the noise scale parameter β provide an interpretable measure of how confidently each constitutive model explains the experimental data. The marginal noise scale distribution, P(β|D,M*), summarizes the effective measurement uncertainty after integrating over plausible parameter combinations. The conditional posterior at the MAP parameter values, P(β|D,M*,θ_MAP), reflects the noise level required when the model is constrained to a single best-fit point. When the red (conditional) curve is noticeably broader or shifted relative to the blue (marginal) curve, it indicates that the MAP parameter set alone absorbs localized mismatches between the model and the data. After marginalizing over nearby parameter values, those mismatches are partially accommodated by variations in θ, producing the narrower blue distributions.

This Bayesian formulation separates variability due to measurement noise from uncertainty in the parameters and from genuine model–data mismatch, making each source of error interpretable. The noise parameter does not compensate for missing physics; it grows only in response to true measurement variability, not model-form error. As a result, the noise scale posteriors provide a transparent measure of confidence in the inferred material response and clearly distinguish between three sources of variability: experimental noise, parameter uncertainty, and model-form error.

When the posteriors coincide, the model performs consistently over the parameter space, whereas when the posteriors are separate, there is a clear parameter set that performs the best. This coinciding behavior does not imply that only one dataset, e.g., UT2 under qKV, is well explained by its selected model. Rather, the contraction or expansion of the noise scale posterior is dataset-specific. For UM1 and UT1, the marginal posterior remains sharply concentrated while the conditional posterior shifts to higher values and is broader. This represents a case where the model can explain the data globally, yet localized misfits at the MAP parameters introduce small, structured residuals.

4.4 Future directions

While the Bayesian IMR method provides a rigorous platform for simultaneous model selection, parameter inference, and quantification of data uncertainty, several extensions remain open for exploration. First, the current implementation relies on strain-based weighting, which prioritizes data points that meet the specified thresholds. Future work could employ adaptive or data-driven weighting strategies, such as variational Bayesian or optimization-based strategies,^96,97 to dynamically balance sensitivity across the different strain-based regimes. Second, although the hierarchical likelihood marginalizes over the noise scale β, the residuals are presently assumed to be Gaussian. Relaxing assumptions through heavy-tailed or mixture models could improve robustness to experimental artifacts and outliers.^88,93 Third, higher-dimensional parameter spaces may benefit from more efficient sampling methods, such as Markov chain Monte Carlo, nested sampling, or surrogate accelerated inference.^69,70 Additionally, in cases where there may not be a clear selected model, Bayesian model averaging can propagate uncertainty in model choice rather than enforcing a single selected model, while cross-validation strategies can assess generalizability across datasets.^98,99 Finally, integrating multi-modal experimental observables, such as digital image correlation fields or bubble asphericity,^37,100 within a joint Bayesian likelihood could strengthen inference in cases where radius–time data alone are ambiguous.

Although the present study adopts the spherical stress integral formulation defined in eqn (2), the Bayesian framework itself does not depend on this reduced representation. Since the technique is independent of the forward solver, it could be extended to consider other constitutive material models. More broadly, this formulation is a transferable method for Bayesian model selection and uncertainty quantification in nonlinear, data-limited systems where understanding rate-dependent constitutive behavior is critical. Thus, the stress integral formulation used in IMR is a case of a broader class of Bayesian model selection problems in computational physics. Other potential applications include the characterization of biological tissues and soft organs under impact or ultrasound stimulation, as well as cavitation- or fracture-induced failure in elastomers and soft composites. Additionally, this method can be used to study high-rate compaction or damage evolution in geological and structural systems. Since the approach integrates physics-informed forward modeling with hierarchical noise treatment and likelihood-based model selection, it can quantify both parametric and structural uncertainty across different parameter spaces.

5 Conclusions

We formulated a hierarchical Bayesian method for inertial microcavitation rheometry of soft hydrogels, combining large-scale forward simulation libraries with statistically grounded model selection. By constructing physically informed priors when evaluating the likelihood across a grid of constitutive models and parameter values, the method discriminates among competing models while providing maximum a posteriori parameter estimates that serve as informed initial guesses. Synthetic data tests show that the method recovers the generating models and parameters and subsets of more complex models can lead to model misidentification if uninformed priors are used. The KV model consistently described the UM datasets and UT1, while UT2 needed the nonlinear stiffening term in the qKV model. The selected models reproduced the collapse and rebound windows, but in some cases, agreement was observed for the first few collapses and subsequently averaged out later noisy oscillations. Posterior predictive intervals and residual diagnostics further confirm that the inferred models reproduce the observed experimental variability across datasets. A comparison between gelatin samples of the same concentration from the University of Michigan and the University of Texas at Austin revealed that the presented method identifies a consistent constitutive model. Models that rely on artificially large noise parameter β values are automatically downweighted in their likelihood, ensuring that plausibility reflects both fit quality and experimental variability. Balancing between simplicity and flexibility, the method offers a ranking of constitutive models, interpretable parameter estimates, and uncertainty diagnostics. By unifying constitutive inference, model selection, and uncertainty quantification, the Bayesian IMR method provides a transferable, probabilistic method for dynamic rheometry of soft matter and polymer materials.

Conflicts of interest

There are no conflicts to declare.

Data availability

The respective Bayesian inference and forward IMR solvers used in this work can be found at https://github.com/vsnchz248/BIMR and https://github.com/InertialMicrocavitationRheometry/IMRv2.git.

Appendices

A Convergence study for grid resolution

We performed a convergence study of the numerical quadrature evaluation of the marginal likelihood. The log likelihood, log Z = log P(D|M_i), was computed using increasing numbers of grid points per parameter axis N, while maintaining the same parameter bounds. For each constitutive model, log Z was evaluated as a function of N using the same likelihood and prior definitions as in the main analysis. Fig. 6 shows the log evidence (top) and relative error (bottom) stabilizing at N = 32, N = 32, and N = 8 for one-, two-, and three-parameter models, respectively. For the relative error calculation, the reference log likelihood value was at the highest resolution for each respective model. The model evidence estimates are converged within 1% tolerance at resolutions of N = 4096^1/d, where d is the model dimension.

Fig. 6 Grid convergence of log model evidence log Z evaluated on synthetic KV data for all seven constitutive models. Top: log Z as a function of grid resolution N per parameter axis. Bottom: Relative error from the finest-grid estimate, |log Z(N) − log Z(N_final)|/max_k|log Z(N_k)|. The dashed line marks the 1% threshold.

B Robustness of the redundancy-based parameter prior

To assess the sensitivity of the redundancy-based parameter prior to reasonable modeling choices, we recomputed the Kelvin–Voigt prior P(θ|KV) under several alternative definitions of the redundancy construction. These tests examined variations in the tolerance scale τ_C, the mismatch norm, and the normalization used in the stress comparison. Specifically, we varied the tolerance scale by multiplicative factors (0.5τ_C, τ_C, and 2τ_C), considered alternative mismatch norms (weighted/unweighted L₂ and L₁), and compared different weighted scalings. In all tested cases, the prior consistently assigns reduced weight to regions where the Kelvin–Voigt model approaches limiting simpler behavior, while retaining higher prior mass where viscous and elastic stresses remain distinguishable. Fig. 7 shows a representative prior map illustrating the qualitative structure of the prior.

Fig. 7 Qualitative structure of the redundancy-based parameter prior P(θ|KV) for the Kelvin–Voigt model. Regions where viscous and elastic stresses become indistinguishable receive reduced prior mass, while regions where both contributions remain distinguishable retain higher prior weight.

To quantify the similarity between these alternatives and the baseline prior used in the manuscript, we computed the Spearman rank correlation between the corresponding prior maps.¹⁰¹ The resulting correlations are summarized in Fig. 8, showing that the prior structure is largely preserved across the tested variants.

Fig. 8 Spearman rank correlation ρ_S between the baseline KV prior used in the manuscript and priors computed using alternative definitions of the redundancy construction. Blue bars correspond to tolerance-scale definitions τ_C, red bars to mismatch norm choices, and yellow bars to weighting schemes. High correlations indicate that the qualitative structure of the redundancy prior is robust to these modeling choices.

Acknowledgements

VS and MRJ acknowledge Dr. Tianyi Chu for the useful conversations in the preparation of this manuscript. MRJ acknowledges support from the U.S. Department of Defense under the DEPSCoR program Award No. FA9550-23-1-0485 (PM Dr. Timothy Bentley). JBE, MRJ, and JY gratefully acknowledge support from the U.S. National Science Foundation (NSF) under Grant No. 2232426, 2232427, and 2232428, respectively. JBE and SHB acknowledge support from the U.S. Department of Defense, the Army Research Office under Grant No. W911NF-23-10324 (PMs Drs. Denise Ford and Robert Martin). This work used Anvil at Purdue University through allocation MCH220010 from the Advanced Cyberinfrastructure Coordination Ecosystem: Services & Support (ACCESS) program, which is supported by U.S. National Science Foundation grants #2138259, #2138286, #2138307, #2137603, and #2138296. Funding agencies were not involved in study design; in the collection, analysis and interpretation of data; in the writing of the report; or in the decision to submit the article for publication. The opinions, findings, and conclusions, or recommendations expressed are those of the authors and do not necessarily reflect the views of the funding agencies.

Notes and references

1. P. A. Taylor, R. Ford and J. S. Ludwigsen, J. Biomech. Eng., 2009, 131, 061007.
2. J. L. Marsh and S. A. Bentil, Front. Neurol., 2021, 12, 626393.
3. E. Yeats, N. Lu, J. R. Sukovich, Z. Xu and T. L. Hall, Ultrasound Med. Biol., 2023, 49, 1182–1193.
4. A. D. Maxwell, C. A. Cain, T. L. Hall, J. B. Fowlkes and Z. Xu, Ultrasound Med. Biol., 2013, 39, 449–465.
5. C. E. Dougan, Z. Song, H. Fu, A. J. Crosby, S. Cai and S. R. Peyton, Biophys. J., 2022, 121, 2721–2729.
6. J. B. Estrada, C. Barajas, D. L. Henann, E. Johnsen and C. Franck, J. Mech. Phys. Solids, 2018, 112, 291–317.
7. J. Yang, H. C. Cramer III and C. Franck, Extreme Mech. Lett., 2020, 39, 100839.
8. S. C. Haskell, N. Lu, G. E. Stocker, Z. Xu and J. R. Sukovich, J. Acoust. Soc. Am., 2023, 153, 237–247.
9. J. B. Estrada, H. C. Cramer III, M. T. Scimone, S. Buyukozturk and C. Franck, Brain Multiphys., 2021, 2, 100034.
10. E. Bar-Kochba, M. T. Scimone, J. B. Estrada and C. Franck, Sci. Rep., 2016, 6, 30550.
11. J. C. Roberts, J. V. O’Connor and E. E. Ward, J. Trauma Acute Care Surg., 2005, 58, 1241–1251.
12. J. D. Finan, T. E. Vogt and Y. Samei, Exp. Fluids, 2024, 65, 114.
13. C. W. Barney, C. E. Dougan, K. R. McLeod, A. Kazemi-Moridani, Y. Zheng, Z. Ye, S. Tiwari, I. Sacligil, R. A. Riggleman and S. Cai, et al., Proc. Natl. Acad. Sci. U. S. A., 2020, 117, 9157–9165.
14. J. Lang, R. Nathan, D. Zhou, X. Zhang, B. Li and Q. Wu, Phys. Fluids, 2021, 33, 031908.
15. B. Rashid, M. Destrade and M. D. Gilchrist, J. Mech. Behav. Biomed. Mater., 2012, 10, 23–38.
16. H. Saraf, K. T. Ramesh, A. M. Lennon, A. C. Merkle and J. C. Roberts, J. Biomech., 2007, 40, 1960–1967.
17. G. Franceschini, D. Bigoni, P. Regitnig and G. A. Holzapfel, J. Mech. Phys. Solids, 2006, 54, 2592–2620.
18. C. R. Siviour and J. L. Jordan, J. Dyn. Behav. Mater., 2016, 2, 15–32.
19. J. Yang and C. Franck, Dynamic Behavior of Materials: Proceedings of the 2019 Annual Conference on Experimental and Applied Mechanics, 2019, vol. 1, pp. 175–179.
20. E. Johnsen and L. Mancia, Journal of Physics: Conference Series, 2015, p. 012022.
21. C. Barajas and E. Johnsen, J. Acoust. Soc. Am., 2017, 141, 908–918.
22. W. W. Chen and B. Song, Split Hopkinson (Kolsky) bar: design, testing and applications, Springer Science & Business Media, 2011.
23. C. van Sligtenhorst, D. S. Cronin and G. W. Brodland, J. Biomech., 2006, 39, 1852–1858.
24. J. Zhang and T. Wierzbicki, J. Mech. Behav. Biomed. Mater., 2011, 4, 1514–1522.
25. M. M. Trexler, A. M. Lennon, A. C. Wickwire, T. P. Harrigan, Q. T. Luong, J. L. Graham, A. J. Maisano, J. C. Roberts and A. C. Merkle, J. Mech. Behav. Biomed. Mater., 2011, 4, 1920–1928.
26. A. Albrecht, K. Liechti and K. Ravi-Chandar, Exp. Mech., 2013, 53, 113–122.
27. H. Khatam, Q. Liu and K. Ravi-Chandar, Acta Mech. Sin., 2014, 30, 125–132.
28. W. W. Chen, J. Dyn. Behav. Mater., 2016, 2, 2–14.
29. B. Song and W. W. Chen, Lat. Am. J. Solids Struct., 2005, 2, 113–152.
30. J. Seifert, K. Sato, M. G. D. Geers, J. A. W. van Dommelen and G. W. M. Peters, J. Biomech. Eng., 2023, 145, 031004.
31. J. Yang, J. Tao and C. Franck, Exp. Mech., 2021, 61, 1181–1191.
32. A. Baldit, D. Ambard, F. Cherblanc and P. Royer, PhotoMechanics 2013, 2013, pp. 4-p.
33. J. Ning, V. G. Braxton, Y. Wang, M. A. Sutton, Y. Wang and S. M. Lessner, Microsc. Microanal., 2011, 17, 81–90.
34. J. Yang and K. Bhattacharya, Exp. Mech., 2021, 61, 719–735.
35. G. Sugerman, J. Yang and M. Rausch, Exp. Mech., 2023, 63, 585–590.
36. B. Pan, K. Qian, H. Xie and A. Asundi, Meas. Sci. Technol., 2009, 20, 062001.
37. A. McGhee, J. Yang, E. Bremer, Z. Xu, H. Cramer III, J. Estrada, D. Henann and C. Franck, Exp. Mech., 2023, 63, 63–78.
38. J. Yang, Y. Yin, A. K. Landauer, S. Buyukozturk, J. Zhang, L. Summey, A. McGhee, M. K. Fu, J. O. Dabiri and C. Franck, SoftwareX, 2022, 19, 101204.
39. S. Buyukozturk, A. Landauer, L. Summey, A. Chukwu, J. Zhang and C. Franck, Exp. Mech., 2022, 62, 1673–1690.
40. P. L. Reu, E. Toussaint, E. Jones, H. A. Bruck, M. Iadicola, R. Balcaen, D. Z. Turner, T. Siebert, P. Lava and M. Simonsen, Exp. Mech., 2018, 58, 1067–1099.
41. P. L. Reu, B. Blaysat, E. Andò, K. Bhattacharya, C. Couture, V. Couty, D. Deb, S. S. Fayad, M. A. Iadicola and S. Jaminion, et al., Exp. Mech., 2022, 62, 639–654.
42. M. Patel, S. E. Leggett, A. K. Landauer, I. Y. Wong and C. Franck, Sci. Rep., 2018, 8, 5581.
43. M. Hrapko, J. Van Dommelen, G. Peters and J. Wismans, Biorheology, 2008, 45, 663–676.
44. G. Puglisi and L. Truskinovsky, Philos. Trans. R. Soc., A, 2016, 374, 20150118.
45. T. Kirchdoerfer and M. Ortiz, Comput. Methods Appl. Mech. Eng., 2016, 304, 81–101.
46. F. E. Bock, R. C. Aydin, C. J. Cyron, N. Huber, S. R. Kalidindi and B. Klusemann, Front. Mater., 2019, 6, 110.
47. K. Upadhyay, J. N. Fuhg, N. Bouklas and K. Ramesh, Comput. Mech., 2024, 1–30.
48. P. Xu, X. Ji, M. Li and W. Lu, npj Comput. Mater., 2023, 9, 42.
49. J. A. Zimberlin, N. Sanabria-DeLong, G. N. Tew and A. J. Crosby, Soft Matter, 2007, 3, 763–767.
50. C. E. Dougan, H. Fu, A. J. Crosby and S. R. Peyton, J. Mech. Behav. Biomed. Mater., 2024, 160, 106698.
51. T. Cohen and D. Durban, Int. J. Eng. Sci., 2010, 48, 52–66.
52. M. T. Warnez and E. Johnsen, Phys. Fluids, 2015, 27, 1–28.
53. R. Gaudron, M. T. Warnez and E. Johnsen, J. Fluid Mech., 2015, 766, 54–75.
54. M. S. Plesset and A. Prosperetti, Annu. Rev. Fluid Mech., 1977, 9, 145–185.
55. C. E. Brennen, Cavitation and Bubble Dynamics, Oxford University Press, 1995.
56. A. K. Abu-Nab, K. G. Mohamed and A. F. Abu-Bakr, J. Phys.: Condens. Matter, 2022, 34, 304005.
57. J. Yang, H. C. Cramer III, E. C. Bremer, S. Buyukozturk, Y. Yin and C. Franck, Extreme Mech. Lett., 2022, 51, 101572.
58. Z. Zhu, S. Remillard, B. A. Abeid, D. Frolkin, S. H. Bryngelson, J. Yang, M. Rodriguez and J. B. Estrada, Soft Matter, 2025, 21, 6717–6734.
59. B. A. Abeid, M. L. Fabiilli, J. B. Estrada and M. Aliabouzar, Ultrason. Sonochem., 2024, 103, 106754.
60. J.-S. Spratt, M. Rodriguez, K. Schmidmayer, S. H. Bryngelson, J. Yang, C. Franck and T. Colonius, J. Mech. Phys. Solids, 2021, 152, 104455.
61. S. Madireddy, B. Sista and K. Vemaganti, Comput. Methods Appl. Mech. Eng., 2015, 291, 102–122.
62. J. T. Oden, E. E. Prudencio and A. Hawkins-Daarud, Math. Models Methods Appl. Sci., 2013, 23, 1309–1338.
63. J. Chiacho, M. Chiacho, A. Saxena, S. Sankararaman, G. Rus and K. Goebel, Int. J. Fatigue, 2015, 70, 361–373.
64. I. Babuška, Z. Sawlan, M. Scavino, B. Szabó and R. Tempone, Comput. Methods Appl. Mech. Eng., 2016, 304, 171–196.
65. M. Muto and J. L. Beck, J. Vib. Control, 2008, 14, 7–34.
66. E. Elmukashfi, G. Marchiori, M. Berni, G. Cassiolas, N. F. Lopomo, H. Rappel, M. Girolami and O. Barrera, Adv. Appl. Mech., 2022, 55, 425–511.
67. D. D. Fitzenz, A. Jalobeanu and S. H. Hickman, J. Geophys. Res.:Solid Earth, 2007, 112, 1–18.
68. B. V. Rosić, A. Kučerová, J. Sykora, O. Pajonk, A. Litvinenko and H. G. Matthies, Eng. Struct., 2013, 50, 179–196.
69. S. Sarkar, D. Kosson, S. Mahadevan, J. Meeussen, H. Van Der Sloot, J. Arnold and K. Brown, Cem. Concr. Res., 2012, 42, 889–902.
70. A. Gelman, J. Carlin, H. Stern, D. Rubin, D. Dunson and A. Vehtari, Bayesian Data Analysis, CRC Press, 2021.
71. J. L. Beck and L. S. Katafygiotis, J. Eng. Mech., 1998, 124, 455–461.
72. L. S. Katafygiotis and J. L. Beck, J. Eng. Mech., 1998, 124, 463–467.
73. H. Rappel, L. A. Beex, J. S. Hale, L. Noels and S. Bordas, Arch. Comput. Methods Eng., 2020, 27, 361–385.
74. T. Most, Identification of the parameters of complex constitutive models: least squares minimization vs. Bayesian updating, in Reliability and optimization of structural systems, ed. D. Straub, CRC Press, New York, 2010, pp. 119–130.
75. T. Chu, J. B. Estrada and S. H. Bryngelson, Comput. Mech., 2025, 76, 431–447.
76. T. Chu, J. B. Estrada and S. H. Bryngelson, arXiv, 2025, preprint, arXiv:2505.13283 DOI:10.48550/arXiv.2505.13283.
77. J. B. Freund and R. H. Ewoldt, J. Rheol., 2015, 59, 667–701.
78. P. P. Nagrani, A. J. Thomas, A. M. Marconnet and I. C. Christov, J. Rheol., 2025, 69, 677–692.
79. A. Rinkens, C. V. Verhoosel and N. O. Jaensson, J. Non-Newtonian Fluid Mech., 2023, 322, 105154.
80. A. M. Sunol, J. V. Roggeveen, M. G. Alhashim, H. S. Bae and M. P. Brenner, arXiv, 2025, preprint, arXiv:2510.24673 DOI:10.48550/arXiv.2510.24673.
81. J. B. Keller and M. Miksis, J. Acoust. Soc. Am., 1980, 68, 628–633.
82. A. Prosperetti and A. Lezzi, J. Fluid Mech., 1986, 168, 457–478.
83. F. Jomni, A. Denat and F. Aitken, J. Appl. Phys., 2009, 105, 0533011.
84. C. T. Wilson, T. L. Hall, E. Johnsen, L. Mancia, M. Rodriguez, J. E. Lundt, T. Colonius, D. L. Henann, C. Franck and Z. Xu, et al., Phys. Rev. E, 2019, 99, 043103.
85. E. Bremer-Sai, J. Yang, A. McGhee and C. Franck, Exp. Mech., 2024, 64, 587–592.
86. InertialMicrocavitationRheometry, inca, GitHub repository, 2025, https://github.com/InertialMicrocavitationRheometry/inca.
87. A. Gandin, Y. Murugesan, V. Torresan, L. Ulliana, A. Citron, P. Contessotto, G. Battilana, T. Panciera, M. Ventre and A. P. Netti, et al., Sci. Rep., 2021, 11, 22668.
88. D. S. Sivia and J. Skilling, Data Analysis: A Bayesian Tutorial, Oxford University, Oxford, 2006.
89. P. Gregory, Bayesian Logical Data Analysis for the Physical Sciences: A Comparative Approach with Mathematica Support, Cambridge, 2005.
90. L. Held and D. Sabanes Bove, Likelihood and Bayesian Inference, Springer Statistics for Biology and Health, 2020.
91. T. Blanchard, T. Lombrozo and S. Nichols, Cognit. Sci., 2018, 42, 1345–1359.
92. E. T. Jaynes, Probability Theory: The Logic of Science, Cambridge University Press, 2003.
93. G. E. P. Box and G. C. Tiao, Bayesian Inference in Statistical Analysis, Wiley Classics Library, 2011.
94. A. Gelman, Bayesian Anal., 2006, 1, 515–533.
95. A. A. Neath and J. E. Cavanaugh, Wiley Interdiscip. Rev.:Comput. Stat., 2012, 4, 199–203.
96. C. M. Bishop, Pattern Recognition and Machine Learning, Springer, New York, 2006.
97. D. J. Rezende, S. Mohamed and D. Wierstra, Proceedings of the 31st International Conference on Machine Learning (ICML), 2014, pp. 1278–1286.
98. J. A. Hoeting, D. Madigan, A. E. Raftery and C. T. Volinsky, Stat. Sci., 1999, 14, 382–401.
99. A. Vehtari, A. Gelman and J. Gabry, Stat. Comput., 2017, 27, 1413–1432.
100. J. Yang, A. McGhee, Z. Tong, G. Radtke, M. Rodriguez Jr and C. Franck, bioRxiv, 2025, preprint DOI:10.1101/2025.02.16.638549.
101. C. Spearman, Am. J. Psychol., 1904, 15, 72–101.

---