The mechanical and sensory signature of plant-based and animal meat

Mechanical testing

We test five plant-based meat products: ham style roast tofurky (Tofurky, Hood River, OR), artesian vegan frankfurter plant sausage (Field Roast, Seattle, WA), signature stadium plant hotdog (Field Roast, Seattle, WA), organic firm tofu (365 by Whole Foods Market, Austin TX), organic extrafirm tofu (House Foods, Garden Grove, CA). For comparison, we also tested three animal meat products: wieners classic animal hotdog (Oscar Meyer, Kraft Heinz Co, Chicago, IL), oven roasted spam turkey (Spam, Hormel Foods Co, Austin, MN), and turkey polska kielbasa animal sausage (Hillshire Farm, New London, WI). Table 1 summarizes the ingredients of all eight products. For each meat type, we test at least n = 5 samples in tension, compression, and shear. Figure 6 documents our sample preparation and our mechanical testing.

Sample preparation

We prepare the samples following our established protocols²⁸. Figure 6b, c, e, f illustrate our tension tests, for which we use a custom 3D printed cutting guide and brain sectioning knife to prepare samples of 1 × 1 × 2 cm³. We super glue the samples to glass microscope slides and wait for 30 min until the glue is fully cured. During curing, we drape a damp paper towel over the samples to keep them hydrated. For the compression and shear tests, we prepare cylindrical samples of 8 mm diameter and 1 cm height using a biopsy punch to extract full-thickness cores from the center of each material, and store the samples in a damp paper tower until testing.

Sample testing

For all three modes, tension, compression, and shear, we test the samples raw and at room temperature at 25 °C²⁸. We perform all uniaxial tension tests using an Instron 5848 (Instron, Canton, MA) with a 100N load cell, see Fig. 6e. We use 3D printed custom grips to rapidly mount and unmount the microscope slides for high throughput testing. Figure 7 shows the part dimensions to create these grips. We mount each sample in the grips, apply a small pre-load, and calibrate the initial gauge length L. We determine the pre-load magnitude for each sample individually. We define pre-load as the minimum load needed to remove any slack, based on visual inspection of a force-displacement curve starting with a fully unloaded specimen, generally on the order of 0.5 N. We then increase the tensile stretch quasi-statically at a rate of (dot{lambda }=0.2 %)/s for t = 50 s, until the sample fails. We perform all uniaxial compression and shear tests using an AR-2000ex torsional rheometer (TA Instruments, New Castle, DE), see Fig. 6d. For the compression tests, we mount the sample, apply a small pre-load, and calibrate the initial gage length L. We determine the pre-load for each sample based on its loading curve, with values on the order of 0.5 N. We then increase the compressive stretch quasi-statically at a rate of (dot{lambda }=0.2 %)/s for t = 50 s, up to a total stretch of λ = 0.9. For the shear tests, we apply a 10% compressive pre-load and calibrate the initial gage length L. We then rotate the upper plate quasi-statically at a shear rate of (dot{gamma }=0.2 %)/s for t = 50 s, up to a total shear of γ = 0.1. To prevent slippage of the samples during the shear tests, we use a sandpaper-covered base plate of 20 mm diameter and a sandpaper-covered top plate of 8 mm diameter.

Front, right, and isometric views of our customized 3D printed device to hold glass slides during tensile testing. Two holders need to be printed to complete the set-up. A standard microscope glass slide fits into the 3 mm slot; it mounts and unmounts easily for high throughout testing. All dimensions are in mm.

Analytical methods and data processing

For each sample and each test mode, we use MATLAB (Mathworks, Natick, MA, USA) to smooth the curves using smoothingspline and SmoothingParam = 1. We interpolate each smooth curve over 21 equidistant points in the ranges (1.0le lambda le {lambda }_{max }) for tension, 1.0 ≥ λ ≥ 0.9 for compression, and 0.0 ≤ γ ≤ 1.0 for shear. For each meat product, we select ({lambda }_{max }) in tension as the maximum stretch in the hyperelastic regime, the loading range within which we observe no visible failure for any of the samples. Finally, we average the interpolated curves to obtain the mean and standard error of the mean for each product and report the data in Table 2.

Stiffness

For each testing mode, we extract the stiffness of each product from the data in Table 2 using linear regression. We convert that the tension and compression data {λ; P₁₁} into strain-stress pairs {ε; σ}, where ε = λ − 1 is the strain and σ = P₁₁ is the stress. We postulate a linear stress-strain relation, σ = E ⋅ ε, and use linear regression to determine the tensile and compressive stiffnesses ({E}_{{{rm{ten}}}},{E}_{{{rm{com}}}}=({{boldsymbol{varepsilon }}}cdot hat{{{boldsymbol{sigma }}}})/({{boldsymbol{varepsilon }}}cdot {{boldsymbol{varepsilon }}})). Similarly, we rewrite the shear data, {γ; P₁₂}, as shear strain-stress pairs {γ; τ}, where γ is the shear strain and τ = P₁₂ is the shear stress. We postulate a linear shear stress-strain relation, τ = μ ⋅ γ, convert the shear modulus μ into the shear stiffness, E_shr = 2 [ 1 + ν ] μ = 3 μ, and use linear regression, to determine the shear stiffness ({E}_{{{rm{shr}}}}=3,({{boldsymbol{gamma }}}cdot hat{{{boldsymbol{tau }}}})/({{boldsymbol{gamma }}}cdot {{boldsymbol{gamma }}})).

Kinematics

We analyze all three testing modes combined using finite deformation continuum mechanics^42,43. During testing, particles X of the undeformed sample map to particles, x = φ(X), of the deformed sample via the deformation map φ. Similarly, line elements of the dX of the undeformed sample map to line elements, dx = F ⋅ dX, of the deformed sample via the deformation gradient, ({{boldsymbol{F}}}={nabla }_{{{boldsymbol{X}}}}{{boldsymbol{varphi }}}={sum }_{i = 1}^{3},{lambda }_{i},{{{boldsymbol{n}}}}_{i}otimes {{{boldsymbol{N}}}}_{i}). Its spectral representation introduces the principal stretches λ_i and the principal directions N_i and n_i in the undeformed and deformed configurations, where F ⋅ N_i = λ_in_i. We assume that all meat samples are isotropic and have three principal invariants, ({I}_{1}={lambda }_{1}^{2}+{lambda }_{2}^{2}+{lambda }_{3}^{2}) and ({I}_{2}={lambda }_{1}^{2}{lambda }_{2}^{2}+{lambda }_{2}^{2}{lambda }_{3}^{2}+{lambda }_{1}^{2}{lambda }_{3}^{2}) and ({I}_{3}={lambda }_{1}^{2},{lambda }_{2}^{2},{lambda }_{3}^{2}={J}^{2}), which are linear, quadratic, and cubic in terms of the principal stretches squared. We also assume that all samples are perfectly incompressible, and their third invariant always remains equal to one, I₃ = 1. The remaining two invariants, I₁ and I₂, depend on the type of experiment.

Constitutive equations

Constitutive equations relate a stress like the Piola or nominal stress P, the force per undeformed area that we measure during our experiments, to a deformation measure like the deformation gradient F²⁸. For a hyperelastic material that satisfies the second law of thermodynamics, we can express the Piola stress, P = ∂ψ(F)/∂F − p F^−t, as the derivative of the Helmholtz free energy function ψ(F) with respect to the deformation gradient F, modified by a pressure term, − p F^−t, to ensure perfect incompressibility. Here, the hydrostatic pressure, (p=-frac{1}{3},{{boldsymbol{P}}}:{{boldsymbol{F}}}), acts as a Lagrange multiplier that that we determine from the boundary conditions of our experiments. Instead of formulating the free energy function directly in terms of the deformation gradient ψ(F), we can express it in terms of the invariants, ψ(I₁, I₂), and obtain the following expression, P = ∂ψ/∂I₁ ⋅ ∂I₁/∂F + ∂ψ/∂I₂ ⋅ ∂I₂/∂F − p F^−t.

Constitutive neural networks

Motivated by these kinematic and constitutive considerations, we reverse-engineer our own constitutive neural network that satisfies the conditions of thermodynamic consistency, material objectivity, material symmetry, incompressibility, constitutive restrictions, and polyconvexity by design^29,44. Yet, instead of building these constraints into the loss function, we hardwire them directly into our network input, output, architecture, and activation functions^45,46 to satisfy the fundamental laws of physics. Special members of this family represent well-known constitutive models, including the neo Hooke³², Blatz Ko³⁰, Mooney Rivlin^33,34, and Demiray³¹ models, for which the network weights gain a clear physical interpretation^44,47. Specifically, our constitutive neural network learns a free energy function that is parameterized in terms of the first and second invariants. It takes the deformation gradient F as input, computes the two invariants, I₁ and I₂, raises them to the first and second powers, (∘) and (∘)², applies either the identity or exponential function, (∘) and exp(∘), and summarizes all eight terms in the strain energy function ψ as the network output. Figure 8 illustrates our network with n = 8 nodes, with the following eight-term free energy function,

where w = [ w₁, w₂, w₃, w₄, w₅, w₆, w₇, w₈ ] and ({{{boldsymbol{w}}}}^{* }=[,{w}_{2}^{* },{w}_{4}^{* },{w}_{6}^{* },{w}_{8}^{* }]) are the network weights. From the derivative of the free energy, we calculate the stress,

where the derivatives of the invariants, ∂I₁/∂F and ∂I₂/∂F, depend on the type or experiment. During training, our network autonomously discovers the best subset of activation functions from 2⁸ − 1 = 255 possible combinations of terms. At the same time, it naturally trains the weights of the less important terms to zero.

Isotropic, perfectly incompressible constitutive artificial neural network with two hidden layers and eight terms. The network takes the deformation gradient F as input and calculates its first and second invariant terms, [I₁ − 3] and [I₂ − 3]. The first layer generates powers of these invariants, (∘)¹ and (∘)², and the second layer applies the identity and the exponential function to these powers, (∘) and exp(∘). The strain energy function ψ(F) is a sum of the resulting eight terms. Its derivative defines the Piola stress, ∂ψ(F)/∂F, whose components, P₁₁ or P₁₂, enter the loss function to minimize the error with respect to the tension, compression, and shear data. By minimizing the loss function, the network trains its weights w and w^* and discovers the best model and parameters to explain the experimental data.

Uniaxial tension and compression

In the tension and compression experiments, we apply a stretch λ = l/L, that we calculate as the ratio between the current and initial sample lengths l and L. We can write the deformation gradient F in matrix representation as

In tension and compression, the first and second invariants and their derivatives are I₁ = λ² + 2/λ and I₂ = 2λ + 1/λ² with ∂_λI₁ = 2 λ − 2/λ² and ∂_λI₂ = 2 − 2/λ³. Using the zero normal stress condition, P₂₂ = P₃₃ = 0, we obtain the explicit expression for the uniaxial stress, P₁₁ = 2 [λ − 1/λ²] ∂ψ/∂I₁ + 2 [1 − 1/λ³] ∂ψ/∂I₂, which we can write explicitly in terms of the network weights w and w^*,

Simple shear

In the shear experiment, we apply a torsion angle ϕ, that translates into the shear stress, γ = r/L ϕ, by multiplying it with the sample radius r and dividing by the initial sample length L. We can write the deformation gradient F in matrix representation as

In shear, the first and second invariants and their derivatives are I₁ = 3 + γ² and I₂ = 3 + γ² with ∂_γI₁ = 2 γ and ∂_γI₂ = 2 γ. We obtain the explicit expression for the shear stress, P₁₂ = 2 [∂ψ/∂I₁ + ∂ψ/∂I₂] γ, which we can write explicitly in terms of the network weights w and w^*,

Loss function

Our constitutive neural network learns the network weights, w and w^*, by minimizing the loss function L that penalizes the mean squared error, the L₂-norm of the difference between model and data, divided by the number of data points in tension, compression, and shear²⁸,

We train the network by minimizing the loss function with the ADAM optimizer, a robust adaptive algorithm for gradient-based first-order optimization using the tension, compression, and shear data for all eight meat products from Table 2.

Best-in-class modeling

Instead of looking for the best possible fit of the models to the experimental data, we seek to discover meaningful constitutive models that are interpretable and generalizable⁴⁸, models that have a sparse parameter vector with a limited number of terms⁴⁹. From combinatorics, we know that our network can discover 2⁸ − 1 = 255 possible models, 8 with a single term, 28 with two, 56 with three, 70 with four, 56 with five, 28 with six, 8 with seven, and 1 with all eight terms. For practical purposes, we focus on the the 8 one-term models and the 28 two-term models, set all other weights to zero, and discover the non-zero weights that characterize the active terms. We summarize the results in a color-coded 8 × 8 error plot, as the average of the mean squared error across the tension, compression, and shear data, and report the parameters of the best-in-class one- and two-term models in Table 1. Motivated by the notable differences in the discovered models for tensile stretches up to 10% versus 35% in Fig. 9, we decide to use the full tension data, not just the 10% stretch.

Discovered one-term models, on the diagonal, and two-term models, off-diagonal, using tensile stretches up to 10% versus up to the peak stress of 35%. All models are made up of eight functional building blocks: linear, exponential linear, quadratic, and exponential quadratic terms of the first and second strain invariants I₁ and I₂. The color code indicates the quality of fit to the tension, compression, and shear data from Table 2, ranging from dark blue, best fit, to dark red, worst fit. The larger stretch range of 35% provides a clearer distinction of the quality of fit for the individual models.

Food texture survey

We prepare bite-sized samples of the eight products, five plant-based, tofurky, plant-based sausage, plant-based hotdog, extrafirm tofu, and firm tofu, and three animal-based, spam turkey, animal sausage, and animal hotdog. We bake all samples in an oven until the plant-based products sre sufficiently warm, and the animal products reach a safe internal temperature for eating. We do not add any sauces or condiments, except for a small bit of oil to prevent excess sticking to the parchment paper. We keep the samples warm until serving. We recruit n = 16 participants to participate in three surveys: the ten-question Food Neophobia Survey³⁵ and the sixteen-question Meat Attachment Questionnaire³⁶, see Fig. 10, and our own Food Texture Survey. We instruct each participant to eat a sample of each meat product and rank its texture features according to our survey. The survey uses a 5-point Likert scale with twelve questions. Each question starts with “this food is …”, followed by one of the following features^21,22: soft, hard, brittle, chewy, gummy, viscous, springy, sticky, fibrous, fatty, moist, and meat-like. The scale ranges from 5 for strongly agree to 1 for strongly disagree. This research was reviewed and approved by the Institutional Review Board at Stanford University under the protocol IRB-75418.

The Food Neophobia Survey uses a 7-point Likert scale with 10 questions so the light blue shaded range goes from 10 (neophilic, open to trying new foods) to 70 (neophobic, not open to trying new foods)³⁵. The Meat Attachment Questionnaire uses a 5-point Likert scale with 16 questions so the light orange range goes from 16 (unattached to eating meat) to 80 (very attached to eating meat)³⁶. The box-and-whisker plots of the minimum, first quartile, median, third quartile, and maximum participant scores are plotted in dark blue and dark orange for the Food Neophobia Survey and Meat Attachment Questionnaire, respectively. The black dots show individual scores.