Characterization of the local mechanical texture of animal meat and meat replacements using multi-point indentation

Over consumption of animal meat is a large contributor to environmental change and should be reduced. This reduction can be achieved by replacing high environmental impact animal meat with an alternative produced from plant-based products or innovative in vitro techniques that have negligible environmental impact. However, the general public will likely consume only meat replacements that mimic the aesthetic qualities of animal meat, with a challenging parameter to tune being the local mechanical response which defines the food texture. In this work, we present a method to characterize the local mechanics of food by using multi-point indentation to spatially measure the local elastic modulus. The resulting heterogeneous mechanical maps are quantified using the established auto-correlation method Moran’s I, which can be used to quantify the resemblance of food texture between samples. The presented technique holds the potential to correlate production parameters such as shear rate and chemical composition with perceived food texture.


Introduction
From 2020 to 2050, the world population will rise by an estimated 2 billion (United Nations and Affairs, 2015). To realize an increase of quality of life and reduce the world's environmental footprint, humanity faces major challenges amongst which is food security. Over the past decades health experts have emphasized the detrimental effects of the over consumption of meat for ethical, health, ecological and environmental reasons (Dekkers, et al, 2018). To decrease meat consumption, meat replacements are being developed with two strategies: either with cultured meat, meat produced by in vitro cell culture of animal cells, or with plant-based products (Bonny et al., 2017). Crucial to increase people's willingness to consume meat replacements, the product must mimic the aesthetic qualities of meat such as size, appearance, taste and texture (Macdiarmid et al., 2016).While mimicking meat flavor is possible, a challenging food quality to tune is texture experienced during mastication which is defined by mechanical properties (Mioche et al., 2003). To date, there are analytical techniques capable of measuring these mechanical properties of food, such as rheology (Fischer et al., 2009), the Warner-Bratzler Shear Force (WBSF) which measures tenderness, the texture profile analysis (TPA) which measures hardness, springiness, cohesiveness, gumminess and chewiness, and tensile testing which measures mechanical anisotropy (Dekkers et al., 2018). However, all of these methods produce results based on bulk measurements averaging over any spatial heterogeneity, requiring local properties to be indirectly inferred.
In this work, we present a methodology that locally measures the mechanical response of meat or meat replacement, an important oral sensing property (Mioche et al., 2003), through the use of multi-point indentation resulting in spatial mechanical 'maps' of elastic modulus. Lastly we analyze the data and introduce a new parameter, the decay length ℓ1

Multi-point indentation and humidity control
The multi-point indentation setup translates in x, y a capacitive load cell using three translation stages and pushes a spherical probe on a sample to locally measure the normal force, as seen in Fig. 1. Further technical instrumental details may be found in a previous publication (Boots et al., 2019). A force-distance curve is obtained from which a local modulus, E x,y , is determined using Hertzian contact theory (Popov, 2010), according to with F the normal force measured by the load cell, R the probe radius, δ the indentation depth and ν the Poisson's ratio of the sample. Using this setup, mechanical maps of several meat types are produced. Mapping a sample placed in air results in poor mapping, since the modulus of the material changes relatively fast due to water evaporation, as samples contain up to 80 vol % of water (Huff-Lonergan and Sosnicki, 2002). Furthermore, the instrument takes approximately 1 h to obtain 50 force-distance curves with an acceptable amount of indentation steps to accurately determine the elastic modulus, thus mapping with reasonable spatial accuracy takes ̃24 h for the samples mapped in this work.
Therefore, to reduce the evaporation rate, meat samples were placed in a custom-built chamber, which contained a piezoelectric water nebulizer, recirculating fan, and dry air source operating at flow rate of 10 L/min , see Fig. 1. The humidity was controlled using a temporal feedback loop from a Thorlabs TSP01 humidity logger. The change of the Young's modulus in time, E x,y (t), for the meat replacement used in this research was found to be δE δt = 10 kPa/h during the measurement when using the humidity chamber set at a relative humidity (RH) of 40%, i.e. a non-humidity-controlled environment. Alternatively, at 90% RH, the δE δt =3 kPa/h and at 95% RH δE δt =1 kPa/h. Controlling the chamber humidity during the measurements decreases the change in sample modulus, but does not completely eliminate water evaporation. Therefore, all further experiments are performed at the maximal achievable 95% RH and room temperature. To correct for this change, the raw data is post-processed, as shown in Fig. 2h-j.
The evaporative change in modulus is determined by measuring the modulus at the four corners of the scanned area before and after the experiment. Since the experiment proceeds row-by-row, data is corrected in accordance with experimental time at which it was measured, shown in Fig. 2i. A second post-process is done for a few outliers in the data by spatially averaging these data points with their immediate neighbours. These outliers are caused by communication errors of the translation stages used in the setup with the computer and are displayed in Fig. 2. The results after this second post-process are shown in Fig. 2a An important factor for mechanical structure is the fiber bundle size. For the meat samples tested in this work, individual fiber bundles can be resolved with the naked eye and are, therefore, on the order of 1 mm. To resolve individual fibers, the used spherical probe radius R = 1 mm and the probe was moved across the surface for a scan area A = 20x15 mm with step size L = 0.5 mm in between measurements. Hence, samples are subsampled as shown in Fig. 2, which was done to limit the change of the modulus due to evaporation. The indentation step size was 20μ m with a maximum indentation depth of 500μ m or when a set force threshold for the meat type was reached by the load cell, the probe was retracted and moved to the next position.

Meat samples
In this research, the following were mechanically mapped: meat replacement normal and orthogonal to a processing direction, tofu, raw beef, raw chicken, cooked beef and cooked chicken. The model meat analogue contained soy protein isolate (SUPRO® 500E IP, Solae, St Louis MO, USA), vital wheat gluten (Roquette, Lestrem, France), NaCl, and demineralized water. The analog has a fibrous and soft texture with a final dry matter content of about 30% wt with dimensions of approximately 30x25x100 mm 3 (w x h x l) and is prepared using the shear cell technology on a couette-type shearing device (Krintiras et al., 2016),Cornet, Edwards, van der Goot and van der Sman (2020). All other samples were purchased at a local supermarket and, unless cooked, measured as acquired after cutting to loadable pieces of sample without any other processing. Cooking samples was done in a pan with sunflower oil, flipping the samples over a few times until they resembled consumer meat.

Results
Finding a length-scale for mechanical heterogeneity from mechanical maps of meat with the multi-point indentation setup, mechanical maps of meat and meat analogues are obtained at the same lateral step size L and total scan area A, shown in Fig. 2. These mechanical maps are

Fig. 1. Experimental setup schematic:
Schematic of the multi-position setup that is used for experiments. The load cell is connected to three translation stages allowing motion in three directions, as indicated by the arrows. The sample is kept in a custombuilt humidity chamber to reduce sample evaporation. Water vapour is generated by a custom built humidity-generator, which has an air inlet and a fan, to ensure constant flow of water vapour. The humidity is monitored continuously with a digital hygrometer. The bottom left inset shows indentation with a probe of radius R, a stress profile for indentation number n = 1 and the step size L of the probe to n = 2. Note that the stress profile is heterogeneous due to the presence of a region of higher modulus within the sample denoted by the outlined region.
used to quantify a heterogeneous length scale, ℓ1 2 , indicative for the local variation in modulus, by using spatial autocorrelation. A powerful example of spatial autocorrelation is Moran's I, a technique widely used to study spatial data analysis in a variety of topics such as agriculture (Donfouet et al., 2017), crime (Chung and Kim, 2019), voting (Robinson and Noriega, 2010), and healthcare (Chen et al., 2020). Moran's I is defined as: where I is the Global Moran's I, N is the number of cells, x is the value of E x,y of the cell of interest, x is the spatial mean of x, w ij is a spatial weights matrix and W is the sum of all w ij . The value of Moran's I scales from 1 when the sample is fully locally correlated to − 1 when fully anticorrelated and 0 when for a random arrangement of the data or decorrelated. The resulting Moran's I for the mechanical maps shown in Fig. 2 are calculated for two w ij with the open source software GeoDa (Anselin et al., 2010), shown in Fig. 2k,l. The spatial weights matrix w ij is crucial in determining how the spatial data is correlated and can be defined in multiple ways. In these results, the maps are analyzed by using a 'Queen' connectivity, such that the side and corner sharing neighbouring cells are included in the analysis. Thus a connected weight with no distance weighting and a neighborhood size of 1 is equivalent to, In addition, the maps are analyzed using w ij of a Queen connectivity combined with an inverse distance weight which for a neighborhood size of 2 is equivalent to,

Discussion
All elastic moduli shown in Fig. 2a-g are calculated from a single indentation using Hertzian contact theory, which implicitly assumes a substrate with a homogeneous elastic modulus. However, real experimental samples are mechanically heterogeneous, included those in this work. An inherent consequence of indentation of a mechanically heterogeneous material is blurring of the normal stress response, which leads to blurring of the mechanical map and, therefore, affects the extracted decay length. This normal stress response is blurred over a region larger than the contact area of the spherical probe with the sample, as is depicted in Fig. 1, and is a function of the probe radius R and indentation depth δ. The size of the indentation area is determined by the contact area of the probe with the surface. In this work, a spherical probe is used, resulting in a contact diameter d c = 2 ̅̅̅̅̅ ̅ Rδ √ , for probes smaller than the thickness of the sample (Chan et al., 2012) as is the case for all data in Fig. 2a-g. For the maximum indentation depth reached 500 μm, the contact diameter d cmax ≈ 1.4 mm; thus, each indentation measurement overlaps with the neighbouring indentation measurement, as d cmax > L and shown in Fig. 1. The choice for step size L and scan area A both affect the total duration of the experiment, as the number of modulus scans n = A L 2 . Given that the sample modulus changes over time, as elucidated in section 2, n should be minimized; A cannot be too big or L too small. A is chosen to comprise several fiber bundles and L is chosen to limit n while still being able to resolve individual fibers, determined by a combination of d c and L. In the case an intrinsic feature size is unknown, R and L can be varied to determine a feature size. In this work, R, L and A are kept constant, such that all results for the decay length can be compared.
The decay length, ℓ1 2 , obtained by Moran's I analysis of the mechanical maps, holds information on the mechanical structure of the sample. In case Global Moran's I decays rapidly, ℓ1 2 decreases and the sample is correlated over a short distance. Hence, a large value of ℓ1 2 is related to a more heterogeneous material and lower value ℓ1 2 is related to a more homogeneous material. This trend is seen in the data in Fig. 2. For example, tofu is a commonly used meat substitute that is known for its homogeneous mechanical structure, and therefore a low decay length is expected. Indeed, Fig. 2m shows that tofu has the smallest ℓ1 2 , implying a higher mechanical homogeneity compared to the other samples. For all other samples, the decay length is dissimilar, indicating a range in the mechanical heterogeneity. For example, ℓ1 2 for beef is smaller than for chicken; this is likely related to the fiber bundle size, which is larger for chicken.
Additionally, meat samples are anisotropic materials with an elastic modulus dependent on the orientation of the sample. Therefore, the cutting plane of the meat affects the measured elastic modulus. In Fig. 2b and c, the sample is cut perpendicular and parallel to the surface, respectively. These mechanical maps clearly show a different structure: In Fig. 2c, fiber patterns are visible, whereas the map in Fig. 2b is dominated by an externally applied temperature gradient in the sample production. The resulting ℓ1 2 captures this heterogeneity with the orthogonal direction being significantly lower than the normal direction, as seen in Fig. 2m.
Lastly, the magnitude of ℓ1 2 is dependent on the choice of the w ij , as shown in Fig. 2m. We chose these two w ij as they are commonly used for spatial correlation analyses, like Moran's I (Anselin et al., 2010); these w ij are isotropic and not directional. For example, the fiber patterns seen in Fig. 2c are thus underrepresented in the current isotropic spatial correlation. By contrast, another choice could be a w ij which contains directional orientation with zeros orthogonal to better correlate only in the direction of mechanical anisotropy in the meat samples thus emphasizes this fiber structure; this is of interest for future investigation.

Conclusion & future outlook
The work shown in this paper aims to provide the meat replacement industry with a methodology to quantify the local mechanical resemblance of animal meat to the produced meat replacement. The presented methodology, using a combination of multi-point indentation and Moran's I autocorrelation analysis, shows that a quantitative decay length, that is indicative of local food structure, can be extracted and directly compared between meat replacements and true meat samples. In future work, the presented elastic moduli maps could be crosscorrelated to, for example, local protein content with a bivariant form of Moran's I to find a relation between local food structure and local material content. This local content information could be obtained using imaging techniques such as Raman microscopy (Yang et al., 2016). In addition, to improve the inherently blurred mechanical maps, samples could be indented with probes of multiple radii, to extract decay lengths that hold mechanical information with reduced or even possibly near-zero blurring when the data is extrapolated to an infinitely small probe size.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.