Small and large deformation rheology on pizza cheese as an example of application to study anisotropic properties of food soft materials

Various rheological techniques were applied to study pizza cheese, as a model anisotropic material. Small Amplitude Oscillatory Shear (SAOS) was combined with Large Amplitude Oscillatory Shear (LAOS) rheology, and the data compared to large deformation tests and confocal microscopy. Pizza cheese was analyzed at various stages of processing: before/during mixing, freshly stretched and after 1 and 2 weeks of storage. Directional sampling and testing were of critical importance for data interpretation. LAOS data were evaluated using a recently developed software that, with visual interactive analysis, including layover plots and similarity maps, allowed to provide a rheological fingerprint of the cheese ’ s process and storage-induced structural changes. This work demonstrated that by using complementary rheological and microstructural techniques it is possible to improve our understanding of the anisotropy at different length scales, and to distinguish the differences in the networks among samples after various processing stages.


Introduction
Many food products rely on a structured anisotropic network to deliver desired textural properties.These foods include stretched curd type cheeses, meat products, as well as novel plant-based alternatives developed to replicate the fibrous networks characteristic of their more traditional counterparts.Anisotropy is used to describe materials that are directionally dependent, meaning that they are characterized by phases with a specific orientation.The anisotropic characteristic can be critical for the functional properties, such as mechanical resistance, stretchability and melting in the case of pizza cheese (Fogaça et al., 2017).The rheological characterization of such anisotropic structure brings additional challenges compared to homogeneous samples, which include the characterization of the anisotropy at various length scales not only in terms of microscopic/mesoscopic structure (Ak & Gunasekaran, 1997;Muliawan & Hatzikiriakos, 2007), but also in terms of texture.
Pizza cheese provides a good model system when developing new approaches to study such materials, since the processing steps leading to the formation of its anisotropic structure are well-known.Pizza cheese can be produced by stretching acidified milk curd whereby visible protein fibers are formed.The acidification can be achieved either through direct addition of acid, the use of a starter culture, or by allowing naturally occurring lactic acid bacteria to acidify the milk (Smith et al., 2018).The stretching process includes kneading and stretching the molten cheese mass either in hot water or brine (Breene et al., 1964) or by direct steam injection (Licitra et al., 2017;Sharma et al., 2016).This process serves to melt the fat and align the protein fibers into a three-dimensional network, creating the final structure responsible for the desirable functional properties of the pizza cheese (Sharma et al., 2016).Processing conditions can influence cheese structures, which are also linked to the type of stretching equipment (Feng et al., 2021;Gonçalves & Cardarelli, 2021), forming and packaging (Bast et al., 2015).It is important to note that due to its structural anisotropy, the mechanical properties of pizza cheese are dependent on the direction of the examination (Ak & Gunasekaran, 1997;Bast et al., 2015;Cervantes et al., 1983).Yet, not all published work reports anisotropy characteristics when evaluating mechanical behavior of stretched cheese structures (Fogaça et al., 2017;Muliawan & Hatzikiriakos, 2007;Olivares et al., 2009).Immediately after stretching, the proteins show large compact structures with large water pockets containing coalesced fat globules.During storage rearrangements occur with water redistributing through the protein network (Dahl et al., 2023).In this work, it was hypothesized that using rheology it is possible to probe the textural anisotropy at various length scales providing a fingerprint of the structural development of pizza cheese during processing and storage.
The structures of Mozzarella, pizza cheese and analogs have been greatly investigated using various rheological and texture techniques (Bähler & Hinrichs, 2013;Banville et al., 2016;Feng et al., 2021Feng et al., , 2022;;Muliawan & Hatzikiriakos, 2007;Sharma et al., 2016), including texture profile analysis and Small Amplitude Oscillatory Shear (SAOS).The latter is very frequently applied to analyze the viscoelastic properties of foods at the molecular length scales.Experiments are carried out within the linear viscoelastic (LVE) region, without introducing disruption of sample's structure.In this regime, the dynamic moduli are independent of the amplitude, yet frequency dependent.However, in many processes, the deformation can be large and rapid, resulting in a nonlinear response from the material.Upon increasing the applied shear during rheological analysis from small to large, the nonlinear response of the material can be investigated and the dynamic moduli become both amplitude and frequency dependent (Muliawan & Hatzikiriakos, 2007;Schreuders et al., 2021).Large Amplitude Oscillatory Shear (LAOS) can be applied to derive additional information on the structural changes of the material, resulting in greater relevance to sensory and technological aspects such as particle size changes, moisture release, cohesiveness, and adhesiveness (Hyun et al., 2011;Wang & Selomulya, 2022).
When a large amplitude sinusoidal strain input is applied, the resulting stress waveform will deviate from a sinusoidal wave, with significant contributions from higher-order harmonics.Up to the third harmonic typically allows for collecting most of the nonlinear information (Hyun et al., 2011;Melito et al., 2013;Schreuders et al., 2021).Higher harmonics, also contain additional information about a material's rheological response, however often do not further contribute to the material characterization.Even harmonics are often associated with experimental errors and therefore usually neglected.However, some studies report the even harmonics to contain additional information such as orientation of the material's structure (Ewoldt et al., 2010;Hyun et al., 2011;Macias-Rodriguez et al., 2018).A three-dimensional space analysis of the dependence of strain, strain rate and stress is then conducted in Lissajous-Bowditch curves (van Der Linden & Sagis, 2001).Within this work these curves will be referred to as Lissajous curves.For stress-strain Lissajous curves, when viewed from the elastic perceptive, a perfect elastic material, storing all energy inputs, will show a straight line.From a viscous perspective, plotting stress vs strain rate, a perfect elastic response appears as a circle.Additionally, from the elastic perspective, a perfect viscous material, dissipating all energy during deformation, will be circular in shape, a plastic material will be characterized as a rectangle and viscoelastic materials will show ellipses.Viscoelastic solids show a short minor axis in the elastic plane and a long minor axis in the viscous plane when viewed from the elastic perspective (Schreuders et al., 2021).
Due to the complexity of the acquired data, LAOS is seldom applied to study the rheological properties of food.However, recent studies have illustrated its potential.For example, it has been used to investigate the structuring of emulsion-filled gels (Ntone et al., 2022), describe droplet-droplet interactions (Fuhrmann et al., 2022), characterize the rheological properties of the plant protein blends as a function of temperature (Schreuders et al., 2021), and in evaluating thermo-mechanical processing of anisotropic plant protein structures (Cornet et al., 2022;Snel et al., 2023).Further applications of LAOS in food networks has been summarized in recent review papers (Joyner, 2021;Wang & Selomulya, 2022).Early studies evaluating nonlinear viscoelastic properties of Mozzarella, pizza cheese and cheddar (Tariq et al., 1998) pointed to LAOS as a useful technique to compare structural features.Sharma et al. (2016) applied LAOS on various commercial cheeses relating their nonlinear viscoelastic properties with sensory and oral processing.Muliawan (2008) used strain relaxation and yield stress measurements to characterize the nonlinear rheological properties of Mozzarella cheese and demonstrated that significant structural changes occur with time, temperature and mode of deformation.More recently, Yildirim-Mavis et al. ( 2022) compared LAOS with tack analysis to understand the deformation-related stretchability of kashar cheeses.Piñeiro-Lago et al. (2023) applied stress-controlled LAOS rheology to acid-curd, to study the thermo-rheological response of the cheese and exploring higher harmonics, local measures, and Lissajous curves in their analysis.Yet, more work is needed to further establish general connections between nonlinear rheological properties and the macroand microstructural features of food materials to advance the interpretation of the physical meaning of LAOS parameters.Furthermore, an improved understanding of how LAOS parameters can be used to characterize anisotropic foods is needed.
Our study objectives include the validation of the LAOS technique by its application to a widely recognized anisotropic system, namely pizza cheese.By utilizing existing knowledge of pizza cheese structure and comparing LAOS results with those obtained through conventional techniques, we seek to enhance our comprehension of the link between LAOS results and structural insights, as well as assess its capability to capture anisotropic characteristics in structures.Samples at three different stages of processing, expected to have distinct rheological properties, were evaluated by LAOS and complementary techniques for characterization of mechanical properties and microstructures; namely confocal laser scanning microscopy (CLSM), texture profile analysis (TPA) and SAOS.In spite of the increasing interest in LAOS rheology, existing software to evaluate LAOS data has not received similar attention, and the current available features can limit the interpretation of nonlinear rheological data, especially when many samples are being analyzed, or structural changes between samples need to be compared.The LAOS results in this work were evaluated using VAOS (Visualization of Amplitude Oscillatory Shear) -a novel open source visualization software specifically developed to better follow changes in structural details in food and soft materials (Madsen et al., 2022).Hence, another novel aspect of this study is the employment of these new visualization tools available in the VAOS software to analyze in detail the LAOS data of pizza cheese at various stages of processing and storage.The improved visual representation of LAOS data is expected to enhance the understanding of the complex LAOS data, by improved comparison of sample responses.

Cheese sampling and preparation
Samples from pizza cheese produced by direct acidification were kindly provided by Arla Foods (Denmark) R&D.Fresh rennet curd was mixed with butter fat and acidified in a pilot scale steam stretcher to reach a final pH of 5.3.Direct acidification was employed for this study to minimize the effect of ripening on structural development.The final cheeses had a composition of 20-25% protein, 14-20% fat, and 45-50% dry matter.Three samples, expected to show differences in their rheological properties were studied: rennet curd (with no fat added yet) (sample 1); a sample of cheese mix in the middle of the kneadingprocessing step, after addition of the fat fraction (sample 2); and a sample of the final cheese immediately after production (sample 3).Furthermore, the changes in the cheese structure with storage were evaluated by analyzing the final cheese (sample 3) after one week (sample 4) and after two weeks of storage at 4 • C (sample 5).All samples were analyzed after cooling, at 10 • C unless otherwise indicated.Samples were cut using a cork borer (20 mm diameter) either longitudinal (L) or transverse (T) to the principal axis of the cheese block as J.F. Dahl et al. illustrated in Fig. 1.This sampling technique has been previously employed in published research (Fogaça et al., 2017;Olivares et al., 2009).For TPA the cylindrical sample sticks were cut to 20 mm in height, 5 mm in height for confocal laser scanning imaging and for rheological analysis slices to thin disks of 2-3 mm in height (Fig. 1).

Texture profile analysis
The texture profile analysis (TPA) was performed on a TMS-Touch texture analyzer (Food Technology Corporation, US) equipped with an ILC-S 1000 N (Mecmesin) load cell and a stainless steel plate probe (7.5 cm diameter).The core temperature of the samples during analysis was kept at 10 • C. The cylindrical samples underwent a double compression test, involving two deformations with a short break in between.The test specifications included a 40% deformation, compressing the samples to 12 mm height to prevent early fracture, a test speed of 200 mm/min, and a trigger force of 0.5 N to ensure accurate measurement initiation.The texture parameters derived from the TPA include peak force, cohesiveness and springiness.Peak force 1 and 2 correspond to the maximum force applied during the 1st and 2nd compression.Cohesiveness defines how a sample withstands the second deformation relative to its resistance under the first deformation, being an indication of resistance to mechanical deformation.Cohesiveness is calculated by the area of work during the second compression divided by the area of work during the first compression.Springiness is defined by the return after the deformation during the first compression, whereafter is has been allowed time to relax between the compressions and can be related to elasticity of samples.Springiness is measured by the distance of the detected height during the second compression divided by the original compression distance (Szczesniak & Hall, 1975).

Rheology
The mechanical response, in both the linear and nonlinear regime, was characterized with a HR20 Discovery Hybrid Rheometer (TA Instruments, US), containing both stress-controlled and strain-controlled modes.The strain-controlled mode provides inertia-free strain-frequency data, avoiding the shortcoming of torque inertia otherwise necessary for LAOS measurements on stress-controlled rheometers (Yazar et al., 2019).Analyses were performed with 20 mm geometry crosshatched parallel plates, the lower being stationary and the upper rotating during measurement.Trimming was applied at a trim gap 10% higher than the test gap to ensure full coverage of the plate surface by the sample.The test gap was controlled by applying a constant axial force of 1N.The temperature of the samples was maintained by the Peltier unit.
Strain amplitude tests were performed at 10 • C from 0.1 to 150% strain (γ) at a constant frequency (ω) = 1 Hz to define the limit of the LVE region and the yield point.The critical strain was defined as the point where the G′ and G″ modules are no longer constant upon increasing strain, corresponding to the end of the LVE region.The critical strain was determined as the intersection between a line fitted to the plateau region and a second line fitted to the drop-off region (Qian & Kawashima, 2016).The materials were exclusively analyzed at a frequency of 1 Hz due to its relevance to the targeted food applications.The yield point, also sometimes referred to as yield stress, is the point on the stress-strain curve that indicates the limit of elastic behavior and the start of plastic deformation, and was determined using the tangent intersection point method described by Willenbacher and Lexis (2018).This method involves drawing two tangents on the graph: one tangent follows the slope of the curve for low strain values, while the second tangent follows the slope of the inflection point beyond the yield point.The inflection point is defined as the point with the minimum slope.The yield point can then be determined by identifying the cross-point of the two tangents.This method is particularly useful for soft or progressively yielding materials where the cross-point of G′ and G″ might be parallel over an extended region.The melting properties were also observed by performing temperature ramps from 5 to 75 • C, with 5 • C/min while keeping constant γ = 1% strain and ω = 1 Hz, being within the LVE region for all samples.The melting point was defined by the cross-over of G′ and G″, when the samples shifted from a solid (G′ > G″) to a fluid (G′′ > G′).
LAOS rheology was carried out with a strain amplitude test, in transient data acquisition mode, with a 3 s conditioning time, resulting in collection of data in 17 separate steps over the strain range from 0.1 to 150%, at constant ω = 1 Hz and 10 • C. The stress responses when exposed to 0.1%, 1%, 10%, 15%, 25%, 40%, 70%, and 150% amplitude were collected and included in the analysis.The cheese after one week of cold storage was further exposed to elevated conditions where amplitude sweeps ranging from 1% to 1000% strain were carried out individually at 5 • C, 10 • C, 20 • C, 40 • C and 80 • C, for evaluating more pronounced nonlinear responses.Data processing and analysis were carried out using the VAOS softwarean open-source visual analysis software (freely available at https://vis-au.github.io/vaos/)offering a specific focus on the comparative analysis across multiple samples and tests (Madsen et al., 2022).The data was converted with Fourier transformation including 300 Points per Quarter Cycle.Using the software, it was possible to inspect the harmonics and evaluate their inclusion.Including the 3rd harmonic was found to have an appropriate signal-to-noise ratio across samples and chosen also to allow comparison to previous work.Axes in Lissajous curve plots reported as unitless strain [− ] ratio versus stress [Pa].Quantitative values of strain-stiffening index (S) and shear-thickening index (T) were derived from the VAOS software, calculated according to equations ( 1) and ( 2), as defined by Ewoldt et al. (2008): G L ' represents the elastic modulus at the largest imposed strain and G M ' the minimum resolvable strain, while η ' L and η ' M represent the apparent viscosity at the largest and minimum resolvable shear rate, respectively.Nonzero values of S and T indicate nonlinearities.A positive S value indicates strain-stiffening properties, a negative value, strain-softening.T captures intracycle viscous nonlinearities, with shearthickening characteristics when positive and shear-thinning when T is negative (Ewoldt et al., 2010;Melito et al., 2013;Yazar et al., 2019).

Confocal Laser Scanning microscopy (CLSM)
Nile red and fluorescein isothiocyanate (FITC) were used as fluorescent staining agents for fats and proteins, respectively.Nile red was applied in 10 μL at a concentration of 0.1 mg/ml dissolved in acetone and FITC 20 μL at a concentration of 1 mg/ml dissolved in acetone.The staining solutions were placed on a microscope slide and acetone was allowed to evaporate before placing the cheese sample and leaving this in contact with the dyes for approximately 10 min before imaging.
Micrographs were obtained by a Nikon C2.5i Eclipse TI2 Racon Confocal laser scanning microscope (Minato, Tokyo, Japan) with a 100 × oil objective in 1024× 1024 pixels sections.The Ar/Kr and HeNe laser beam were set at 488 nm and 561 nm, respectively.

Statistical analyses
Statistical analyses were carried out in R-studio (2021.09.2).Oneway analysis of variance (ANOVA) and TukeyHSD analysis were conducted to test for significant differences between TPA parameters and the quantitative rheological measures (critical strain, yield point and melting point) and LAOS (stiffening and thickening index and maximum stress).TukeyHSD analysis was performed to test for significant difference between groups.

Microstructure
Fig. 2 shows the microstructure by CLSM of the various cheese samples sampled either longitudinally or transversely to the protein fibers.The curd (Fig. 2a) showed an isotropic dense protein network, as fat was not yet added at this stage.After addition of fat and heating with steam, at the stage of kneading and stretching (Fig. 2b) the cheese showed larger strands oriented in the direction of the stretching.At this stage, the fat still shows spherical domains in large, hydrated areas, which in the final cheese (Fig. 2c-e) are much reduced.During the storage period of one (Fig. 2d and h) and two (Fig. 2e and i) weeks, diffusion of fat domains and swelling of the protein network occurs.Structural anisotropy was clearly observed for samples cut longitudinally to the protein fibers (Fig. 2b-e), with larger and more elongated fat regions interlocked in the protein network.In contrast, transverse cut samples showed little anisotropy (Fig. 2f-i), with smaller spherical slices elongated fat domains.This observation emphasizes the significance of considering the sampling direction in the analysis of anisotropic structures, as it is likely to impact the material's mechanical response during large deformation studies.The microstructural observations shown in Fig. 2 are in agreement with previous literature reports (Auty et al., 2001;Bast et al., 2015;Feng et al., 2021;Ma et al., 2013).

Texture profile analysis
Texture profile analysis (TPA) is often used to evaluate the large deformation properties of foods as it relates to some important texture sensory attributes (Foegeding et al., 2011).Average calculated values of the TPA parameters are included in Fig. 3 for samples cut either longitudinally (L) and transversely (T) to the internal protein fibers.Significant differences between samples were observed for peak force at the first compression (p < 0.05), with the maximum force increasing from the molten mix to the fresh cheese and then stored for 1 week.This clearly showed structure development during processing of pizza cheese.No significant difference is detected between 1 and 2 weeks of storage, indicating that after the first rearrangements, hardness is less affected by the changes during cold storage.Cohesiveness and springiness were found to be the highest for the curd (sample 1, without fat), due to its very dense protein matrix, yet not found significantly different.No significant differences were observed in the TPA parameters between sampling direction (p < 0.05) apart from sample 2 and 4, where an increased peak force was detected for cheese longitudinally sampled compared to those cut transversely.For sample 2, during the kneading and stretching process, the highest cohesiveness values were observed in longitudinally cut samples, but not found significant.Yet, this finding suggests that the structure is more resistant to deformation when subjected to compression parallel to the protein fibers, but also indicates that sampling direction is a challenge.No significant differences were seen between samples for cohesiveness and springiness.Previous authors (Yildirim-Mavis et al., 2022) have also shown no significant changes in springiness between different cheeses.While peak force measurements are directly related to sample hardness, cohesiveness and springiness, are empirical parameters which depend on more complex properties, and not only dependent on microstructural changes.

Rheological properties
Strain amplitude sweeps were performed to evaluate rheological properties at 10 • C. Fig. 4 shows an overview of the stress vs strain at the end of the LVE region, estimated by the critical strain.Such overviews are of interest as they can distinguish the properties of a material, from brittle to rubbery, and from soft to tough (Truong & Daubert, 2001;Tunick & Van Hekken, 2002).Overall, the curd (sample 1) appears as a strong material with a tough texture, observed by a higher stress response and higher critical strain, linked to the dense protein matrix with no fat.With addition of fat (sample 2) and heating with steam, the sample becomes softer, observed by a lower stress response at the critical strain.Once the network of the cheese is fully created the structure is more prone to deformation, but also more resilient, rubbery, due to the formation of the elastic network.
Analysis of the yield point, in terms of the limit of the elastic range and start of viscous properties, showed that the elasticity range increases from a yield point of 21% strain for the curd sample, to 26% strain for the sample during stretching and 46% strain for the final cheese structure, remaining stable throughout storage (sample 3-5).Both yield point and critical strain point were independent of sampling directions (longitudinal and transverse) indicating that these oscillatory rheological parameters do not capture differences in fiber direction.
In addition to viscoelastic measurements carried out at 10 • C, temperature ramps were performed to evaluate differences in the melting properties.The curd and kneaded mixture exhibited a melting point of approximately 50 • C, which was lower than that observed in the final cheese samples, with a melting point range of 69-71 • C.This variance  can be attributed to compositional differences, the dynamics of the mixing and the equilibration occurring in the stretcher after the addition of the acid to the final pH of 5.3.Furthermore, distinctions in melting points can be ascribed to variations in water mobility and protein alignment, which is not fully developed for the mixture during kneading (sample 2).Heating and stretching the curd increases the proportion of mobile water (Cais-Sokolińska et al., 2018).Hence, more mobile water (in the kneaded mix) as also visible from the CSLM images (Fig. 2), can lead to a lower melting point as the water molecules can facilitate structural breakdown in the transition to a liquid state.No significant differences were found for the melting points among the final cheese samples during storage.These melting points were slightly higher than previously reported melting points for low moisture Mozzarella type cheese, found to range from 50 to 65 • C (Muliawan & Hatzikiriakos, 2007;Sharma et al., 2016), possible due to differences in composition.
In oscillatory rheology usually only information derived from the first harmonic is included.In LAOS rheology higher order harmonics are included in the analyses during in situ dynamic sample deformations.Fig. 5 shows the contribution from the harmonics expressed in normalized power.At 0.1% strain (Fig. 5a), the response is within the linear regime, with clear dominating contribution from the 1st harmonic.With increased strain, however, larger contributions are present from higher odd harmonics, as seen for 100% strain (Fig. 5b).The contribution of the various harmonics was investigated using a heat map as provided by the VAOS software (Fig. 5c).For all samples used in this study the largest contribution is seen for the 1st harmonic followed by the 3rd harmonic, in agreement with previous reported literature (Melito et al., 2013;Hyun et al., 2011).A small contribution from the 5th harmonic could be argued, whereas higher harmonics are seen not to contain additional information.This is in agreement with previous literature, concluding that for cheese, the highest significant harmonic does not normally exceed five (Tariq et al., 1998).The non-linear contributions coming from the instrument itself are typically around 10 − 3 of the fundamental harmonic intensity (Minale et al., 2022).The instrument sensitivity is indicated by the dotted line in Fig. 5a and b.During the experimental setup of this study, oscillatory measurements were initially attempted using smooth parallel plates, which resulted in the cheese samples slipping and an observed dominance of even harmonics in the signal.However, a change in the geometry of the plates to cross-hatched parallel plates led to the disappearance of even harmonics.The absence of contributions from even harmonics, as evidenced by the normalized power contribution being negative for all even harmonics in Fig. 5, suggests that this modified experimental setup had limited slip effects.Other studies have indicated that microstructural information can also be derived from the even harmonics (Hyun et al., 2011), but this does not seem to be the case here.However, further  investigation of even harmonics contribution to degree of anisotropy is of relevance.
Fig. 6 depicts the repeatability of the LAOS approach with Lissajous curves for triplicate measurements.Larger variations between triplicates are observed in the nonlinear regime (150% strain) compared to the ellipse formed in the linear regime (1% strain).Data variability for LAOS analysis of cheese has previously been reported due to the heterogeneous structures and the coalescence of fat globules causing wall slip in the measurements (Tariq et al., 1998).However, the variation in repeatability observed in this work is not associated with wall slip, as a top-and bottom crosshatched geometry was applied and there was minimal contribution from even harmonics (Fig. 5).The decreased repeatability seen in the nonlinear regime was caused by the large deformation occurring differently for individual measurement, and differently depending on the sample.The largest variations were noted for the curd and the mix at kneading (samples 1-2) compared to those of final cheeses which were less heterogenous (samples 3-5), and where the network was further developed.Only the Lissajous curves from the elastic perspective were included in the analysis due to the solid state of the samples.Evaluation of the viscous properties is suggested to be more relevant when analyzing samples in a molten state, as done in recent studies (Yildirim-Mavis et al. 2022).
To follow the evolution of the mechanical response across the strain amplitude, a representative Pipkin diagram was constructed.Fig. 7 illustrates how the shape of the normalized elastic Lissajous curves changes with increasing strain amplitude for the five different cheese samples prepared longitudinally.Overall, an ellipse-shaped curve, representing a viscoelastic response, dominates the diagram.These results agree with previous reports showing similar elliptically shaped loops for Mozzarella at a strain below 50% strain (Sharma et al., 2016).The area enclosed by the curve increased at higher strain values, due to sample deformation.A trend, which is supported by previous studies (Yildirim-Mavis et al. 2022;Tariq et al. 1998) suggesting that the elastic properties progressively changed to more viscous properties, as noted by the wider ellipsoidal shape.This was especially observed for sample 2. The Lissajous curve area is related to energy dissipation (Ewoldt et al., 2010).During the stretching step (sample 2), water and fat are yet to be fully mixed.Fat domains appear more globular based on CLSM images Fig. 6.Repeatability of Lissajous curves from the elastic perspective at 1% and 150% strain amplitude for LAOS measurements for the curd, the sample during kneading and stretching, fresh cheese, cheese stored 1 week and cheese stored 2 weeks.(A-J) Longitudinal (K-R) Transverse sampling.Colors represent three individual performed measurements.Notice differences in scales.(Fig. 2), the cheese matrix demonstrates lower viscosity and greater deformability, leading to a larger Lissajous curve area.The subsequent reduction in area for the final cheese samples indicates a transformation from a more deformable state to a more structured.Upon final processing, where the kneading has been completed, the cheese is a more elastic system that can better withstand deformation.Similar trends in responses were observed for transverse samples, yet with less distinct differences (data not shown).If frequency-dependent data is considered, the use of Pipkin diagrams to depict other dimensionless groups such as the Deborah and Weissenberg numbers, could be relevant in further characterization.
However, one of the main limitations to the use of Pipkin diagrams is the necessity to normalize data for the curves to be at a reasonable size for exploration, making Lissajous responses look much alike across samples and treatments.Plotting overlaying Lissajous curves without normalization enables a more detailed examination of deformation dynamics, facilitating the identification of variations between samples.The VAOS software offers a convenient feature for creating 2D overlay plots of Lissajous curves (see Fig. 8) directly within the software, easing the data analysis process.At low strain values (Fig. 8a), when still within the linear region, samples appeared similar to one another, with narrow ellipse-shaped curves as also evident from the Pipkin diagram.Upon increasing the strain, the area enclosed by the curve increased, reflecting the nonlinear response.Using the 2D overlay plots, the differences in rheological properties between samples during deformation become evident.For example, in the nonlinear regime (Fig. 8b), the curd sample (L1) showed a significantly different behavior compared to that of the remaining samples, due to different composition of the matrix.The sample during kneading (L2), was also different from the final cheese samples.The 2D layover plots make these differences appear more distinct than the Pipkin diagram.There were limited changes during storage (samples L3-L5), yet some differences in the rotation of the curves can be identified, pointing towards slight structural differences.No pronounced differences were observed upon comparing longitudinal samples and transverse (data not shown).However, 2D layover plots for samples probed transversely to the protein strands where less distinct from one another.This indicated that structures can be more uniquely characterized when tested in the longitudinal direction to the protein strands, where the oscillatory shear elongates and bends the fibers, as illustrated in Fig. 1.It is important to mention that the lack of difference between sampling direction might be linked to the small sample size prepared for the rheological measurements.Depending on the fiber length, fibers might slide against one another during deformation, making the length scale of the measurement, and the ability to capture differences, a function of the fiber length.
Following the evolution of the elastic Lissajous curves, the ellipse-formed shape is retained upon increasing strain, as discussed from the Pipkin diagram (Fig. 7).However, comparing the 2D layover plot at small deformation (Fig. 8a) to at large deformation (Fig. 8b) it appears that the ellipse rotates upon the increased strain.A clockwise rotation of the curves upon increasing stain amplitude indicates gradual softening (Ewoldt et al., 2007), which is observed for sample containing fat (sample L2-5).Additionally, a distortion from the elliptical shape indicate strain-stiffening behavior (Ewoldt et al., 2007), which is observed to a larger extent for the sample during kneading (L2).A recent study by Piñeiro-Lago et al. (2023) introduced the concept of Pseudo-Linear LAOS behavior, where the Lissajous curves rotate as amplitude increases while maintaining a low level of distortion.In all types of cheese tested in that study, a pseudo-linear LAOS behavior was observed where rotation was strong and distortion was minimal, indicating that the first harmonics changed significantly in a regime where higher harmonics were small enough for the Lissajous curves to be nearly elliptical.These findings are in agreement with observations of this work for the final cheese structures (sample L3-5).Despite the cheese becoming softer with increasing stress amplitude, the intracycle stress response, being the elliptical Lissajous curve, was maintained with increasing strain.Additionally, Piñeiro-Lago et al. ( 2023) introduced a quantitative metrics, including rotation and distortion, providing a benchmark for differentiating between different types of materials, which will be of interest to explore in future studies.Quantification of the mechanical responses is required to identify how different responses are.
To quantify the evolution of the elastic Lissajous curves in this study, the dimensionless strain-stiffening index and shear-thickening index of nonlinearity are calculated (Equations ( 1) and ( 2)).As shown in Fig. 9, nonlinearities appear at about 40% strain, becoming dominating at 100% strain, in full agreement with the yield point.Longitudinal cheese samples L3-L5 (final cheeses, at varying storage times) showed strainstiffening behavior (Fig. 9a), with increased expansion (elongation) of the protein networks, in accordance with previous reports (Melito et al., 2013;Yildirim-Mavis et al., 2022).The cheese mix (sample L2), collected during kneading and stretching, showed a larger stiffening index (S) than the remaining samples at both 40% and 100% strain amplitude, indicating high elastic nonlinearities for this sample.The larger strain stiffening behavior for sample L2 indicates that the cheese network becomes more strain stiffening as it undergoes structural transformations, as the degree of alignment is strongest during the stretching processing.This more distinct degree of alignment for sample L2 is also visible in the confocal images (Fig. 2).Additionally, the moisture content is expected to play a crucial role in its mechanical behavior, as also discussed during the interpretation of the melting points results.In the final cheese, moisture is expected to be more uniformly distributed and bound within the cheese matrix.This binding of Fig. 8. 2D Layover plot of Lissajous curves from the elastic perspective for pizza cheese during processing and storage at 1% (a), 150% (b) strain amplitude.Cheese samples were cut in the longitudinal direction to protein fibers.Sample key: L1 = curd (dark blue), L2 = during processing (purple), L3 = fresh cheese (green), L4 = cheese stored 1 week (brown) and L5 = cheese stored 2 weeks (light blue).moisture can contribute to a more stable and less deformable structure, reducing strain stiffening of sample L3-L5 compared to sample L2.The curd sample (L1) had values for both S (Fig. 9a) and thickening index (T) (Fig. 9b) closest to zero, demonstrating less nonlinearities compared to the other treatments, linked to no anisotropy structure present in the curd.Anisotropic materials which are more resistant to deformation along their aligned axes may exhibit more pronounced strain stiffening behavior compared to an isotropic system as the curd.The negative values of T (Fig. 9b), observed for most samples, are associated with shear-thinning (pseudoplastic) responses (Ewoldt et al., 2010).Shear-thinning responses were likewise reported by Yildirim-Mavis et al. (2022).The S and T indexes calculated at 40% and 100% strain amplitude were different from the response at the lower strain values (0.1, 1 and 10% strain), for all samples (L1-L5).No significant differences for these stiffening or thickening indexes were found between samples cut in different directions (L vs T).
Further exploration of quantitative measures makes it possible not only to observe pronounced differences between samples, but also more subtle changes.When plotting the maximum stress measured as a function of strain (Fig. 10), samples taken during processing (1 and 2) appear different with higher maximum stress values, and stiffer than the final cheese (Fig. 10a), in agreement with the traditional oscillatory rheology data.In addition, plotting the maximum stress as a function of strain, differences between cheese samples during storage may also become evident, yet not found significantly different.For example, a lower maximum stress is detected for sample 3, analyzed at the day of production, at 70-100% strain, compared to sample 4 and 5, both when analyzing longitudinal (L) and transverse (T) to the protein fibers.The maximum stress quantifies the rotational differences observed in the Lissajous curves at 150% strain (Fig. 8b).In these samples, no proteolytic activity developed because of direct acidification of the cheese, and therefore the change could be attributed to the hydration of the protein and diffusion of the fat, as shown in the microstructure (see Fig. 2).Several studies have investigated the hydration of the protein phase in pizza cheese during storage, by determination of expressible serum (Auty et al., 2001;Guo & Kindstedt, 1995), time-domain NMR T2 relaxometry (Vermeir et al., 2019) or with confocal Raman microscopy (Dahl et al., 2023).
The cheese sample stored for one week was further evaluated at elevated conditions, namely, higher strain values (100-1000% strain) and increased temperatures (10-80 • C), where more variations in rheological responses are expected.No contribution from even harmonics were observed even at these high strains and temperatures, validating that no slip occurred during measurements.2D Area layover plot (Fig. 11a and b), a novel feature of the VAOS software, constructed by overlaying half of the Lissajous curves, refining the visual components, and incorporating a contemporary aesthetic.Unique deformation responses are detected at these extreme testing conditions with smaller maximum stress values in the combination of high strain and high temperature.From the Pipkin diagram (Fig. 11c), it can be derived that at 100% strain and low temperatures (<20 • C), the viscoelastic response still shows an ellipse-shaped curve.Upon increasing temperature (>20 • C) more nonlinear responses, and almost plastic responses are observed once approaching and passing the melting point (>60 • C).Reaching high strain values (500-1000%) also induced a nonlinear response even for samples tested at low temperature.Hence, following the dynamics at elevated conditions may provide additional information in the evaluation of processing-induced changes in complex matrices such as cheese.
An additional novel feature of the VAOS software includes the opportunity to compare responses in a similarity map (Fig. 11d).To this end, pairwise similarities between all elastic Lissajous curves are computed using Dynamic Time Warping (Müller, 2007).Then, one symbol per curve is positioned on the similarity map using force-directed placement (Kobourov, 2014), so that their pairwise distance is inversely proportional to their computed similaritythe closer two symbols are, the more similar are the shapes of the curves they represent.As a consequence of this placement method, absolute positions of the symbols are not meaningful and hence no axes are shown.To better aid in judging relative positions (i.e., distances between symbols), a perfect elastic response serves as a reference point.The closer the responses are positioned to the perfect elastic response (Fig. 11d, red dot), the closer they are to being able to recover their structure.This novel diagramming approach based on the full shape of the Lissajous curves has proven to be more robust than comparing samples based on single point observations, while being highly scalable to showing hundreds of curves at a single glance.Overall, the similarity map highlights the power of performing LAOS testing.Fig. 11d compares curves for one pizza cheese sample, measured at different strains and temperatures.The material response when exposed to low strain values (100%), and different temperatures, is similar.However, at higher strain values (e.g., 500%) and different temperatures suggest a less elastic response.The opportunity of testing at elevated strain values shows the potential to separate the responses.This is clearly recognized by looking at the similarity map for the cheese analyzed at 5 • C. The responses are distributed in a well-defined line pointing away from the perfect elastic response, with increasing strain.Hence, below the melting point, responses follow the expected trend with positions further away from the perfect elastic response when exposed to higher strain values.Above the melting point (≥60 • C), responses appear more clustered, with a less defined trend.The decrease in elastic behavior at these elevated temperatures is due to softening and mobilization of the protein matrix, because of reduced protein-protein interactions and melting of the lipid domains.Yildirim-Mavis et al. ( 2022) have also elucidated the alterations in the microstructure of cheese upon exposure to elevated temperatures.They further explicate that the fat contribution to the elastic modulus decreases faster than the protein matrix contribution at increased temperatures.The similarity map has the ability to capture the compelling differences in sample structure in response to temperature and strain.The evaluation of a particular point in a process, linked to changes in the properties of the product, can be investigated in great detail using this diagramming technique.The results shown in Fig. 11 demonstrate the potential of the new visualization opportunities utilized in VAOS to characterize the details of materials during LAOS deformations and derive important information to design materials of high functionality.

Relating LAOS to complementary techniques
The study of large deformation of soft matter can be approached using both TPA and LAOS rheology.TPA is limited and provides empirical parameters.In contrast, oscillatory rheology is a dynamic test that measures the response of a material to a sinusoidal deformation, offering a wider range of deformation conditions, including precise control of temperature.However, each technique provides valuable information at different length scales, as indicated in Fig. 1.TPA offers insights at the macroscale, allowing the visualization of anisotropy.Hence, anisotropy in pizza cheese correlates with its texture and mouthfeel.CLSM aids in understanding the structural aspects at the microscale.Oscillatory rheology, on the other hand, provides details on viscoelastic properties at smaller length scales.Comparing the mechanical responses derived from TPA and oscillatory rheology illustrates Fig. 11.LAOS rheology at higher strain levels and higher temperature conditions.2D Area Layover plot (from the elastic perspective) for pizza cheese cold stored 1 week at 200% (a) and 1000% (b) amplitude.Fingerprint Lissajous curves (from elastic perspective) (c) and similarity map (d) tested at temperatures ranging from 5 to 80 • C for 1-1000% strain.The similarity map places one symbol per Lissajous curve according to their pairwise distancethe closer two symbols, the more similar they are and vice versaforming clusters of similarly shaped curves and thus of similar material properties.
the complementary nature of the techniques, concluding that length scales of fiber can affect the results of the measurements.Sampling needs to be well-controlled, as directional sampling appears as a challenge for both TPA analysis and oscillatory rheology measurements.The capability of LAOS rheology to capture anisotropy and recognize variances between the movement of fiber planes when subjected to transverse sampling as opposed to bending them during longitudinal sampling (as illustrated in Fig. 1) is anticipated to be directly correlated with fiber size, being a limiting factor.However, the unique capability offered by LAOS to assess structural characteristics under relevant process conditions enables a more comprehensive tracking of deformations taking place during the stretching and kneading process operations.Insights of strain-stiffening and maximum stress behaviors are argued to correlate with degree of anisotropy and water mobility.The new insights provided by LAOS rheology include the opportunity to follow the structure at relevant process conditions, understanding the development of the elastic structure and obtain a structural understanding of large deformation based on mechanical behavior at microscale.
The starting point of the cheese production, the curd (sample 1) made by a concentrated protein fractions before addition of fat, was found unable to form a homogeneous mass, resulting in low peak force values (TPA), a low yield point, low stiffness (LAOS) and a low melting point.Yet, the concentrated protein fractions also showed to be a tough network, withstanding deformation, evidenced by a high cohesiveness value (TPA), high critical strain value and continuous ellipse form responses with increasing strain (LAOS).
Once fat and water are introduced with kneading, and the curd is heated and stretched (sample 2), the matrix becomes softer when exposed to deformation.Yet, increased strength and stiffness are observed during LAOS.Samples taken during the stretching step, and measured transverse to the protein strands, showed the lowest peak force and cohesiveness from the TPA results, indicating the influence of sampling direction on measurement responses.These samples still showed large water domains containing distinct lipid areas in CSLM images (Fig. 2).
Only once the protein network was fully formed and stretched in the fresh cheese (sample 3) the sample elasticity was extended, as seen both from the higher yield point and the ellipse-shaped curve that is less deformed compared to sample 2 at high strain values.This may be caused by a higher resistance of the protein domains when fully developed.Indeed, it was not possible, by evaluating the springiness parameter of the TPA, to identify significant differences among the cheese samples during processing.The oscillatory rheological evaluation of the cheese structures during storage showed the structure to be softer (towards "mushy" in the texture map), in contrast to the TPA high peak force response.This discrepancy illustrates how sample scale can affect the measurement: compressing the protein strands as in the TPA measurements provides a different mechanical response than bending the fibers or rotating the fiber planes during oscillatory measurements.
With LAOS measurements it was possible to capture even subtle changes in complex food networks during processing, highlighting the potential to be used in process control and optimization, as the mechanical response within the nonlinear regime can be explored in great detail.Indeed, although traditional techniques (TPA, SAOS) showed similar properties for the final cheese as a function of storage, subtle textural differences can be better noted using LAOS.Further information can be derived from exploring the LAOS parameter of maximum stress, quantifying differences in the Lissajous curves.The analysis during storage, for example, may suggest that careful determination of maximum stress may highlight microstructural changes linked to an expected hydration of the protein phases.It is therefore possible to conclude that LAOS is a complementary technique to the wellestablished traditional approaches to evaluate the mechanical behavior of pizza cheese, and more in general of foods and soft materials.With the successful validation of the LAOS methodology on a wellestablished anisotropic structure, it can be effectively applied in the analysis of unfamiliar structures, such as plant-based analogs, to gain valuable insights into their mechanical properties and behavior.

Conclusions
Large Amplitude Oscillatory Shear (LAOS) rheology was shown to be able to probe the structure of samples of pizza cheese, an example of an anisotropic material.Using a newly developed software, VAOS, LAOS responses of different cheeses were examined.This allowed us to explore differences in behavior due to stretching and structural changes during storage, by analyzing higher harmonics, local measures, and Lissajous curves, resulting in more advanced structural interpretations compared to traditional mechanical analysis.
Samples of the pizza cheese taken during processing were found significantly different from the final cheese structure in their nonlinear response, indicating that LAOS introduces new opportunities to characterize structured materials as they are being formed.Especially differences in stiffening behavior could be related to degree of anisotropy and water mobility during the stretching step.The pizza cheese samples at different storage times showed similar mechanical responses, indicating a more stable structure compared to the samples collected during processing.Yet, some differences in the viscoelastic properties upon storage were revealed by LAOS and attributed to the microstructural rearrangements.
It is important to point out that the direction of sampling is critical when testing anisotropic structures, and great care needs to be taken not only to prepare the samples for analysis, but also in the data interpretation, based on the type of deformation applied.The inclusion of complementary methodological techniques has shown improved understanding of anisotropy at different length scales, due to differences in investigation direction, i.e., texture profile analysis (TPA) performing large deformation by compressing fiber planes, while oscillatory rheology rotates the fiber planes during deformation.Though, the length of the strands may be a critical factor, depending on the sample length scale, leaving quantitative responses a challenge.Further investigation is suggested to include of frequency dependence data, which might reveal a greater understanding of the network, as it was recently found that anisotropic systems with different energy interactions exhibit different frequency dependencies in their rheological response, whereas isotropic systems exhibit identical energy interactions (Iyer, 2022).

Fig. 1 .
Fig. 1.Illustration of cheese sampling and testing directions for texture profile analysis (TPA), confocal laser scanning imaging (CLSM) and oscillatory rheology.Samples were cut in cylindrical tubes either longitudinal (L), being parallel, or transverse (T), being perpendicular to the principal axis of the cheese block.Below each technique is a schematic illustration indicating the corresponding length scale.From left to right, the length scales are macroscale (yellow and blue areas indicating phase separation), microscale (protein network interlocking elongated fat globules in red), and nanoscale (aligned protein particles in gray supported by calcium ions in black).

Fig. 2 .
Fig. 2. Microstructural images of samples from pizza cheese production obtained by CLSM for sampling longitudinal or transverse to the protein strands and analyzed in the direction of perpendicular cuts: Curd (a), during stretching (b,f), cheese analyzed at the day of production (c,g), cheese stored one week (d,h) and cheese stored 2 weeks at 4 • C (e,i).Green corresponds to protein and red to fat.

Fig. 3 .
Fig. 3. Texture profile analysis parameters: Peak force of first (a) and second (b) compression, Springiness (c) and Cohesiveness (d).Pizza cheese sampled longitudinal (L, solid bars), or transverse (T, open bars), to protein fibers.Standard deviation based on triplicate samples and superscripts indicates significant differences between samples.

Fig. 4 .
Fig. 4. Texture map plotting stress versus strain at the critical strain, defining the end of the LVE region, for pizza cheese during processing and storage, sampled longitudinally.Square = curd, triangle = during stretching, closed circle = cheese at the production day, open circle = cheese stored 1 week and half-filled circle = cheese stored 2 weeks.Values are averages of triplicate samples, with bars indicating standard deviations.

Fig. 5 .
Fig. 5. Contribution from harmonics in normalized power of the contribution intensity for (a) 0.1% strain amplitude, (b) 100% strain amplitude and in a (c) heatmap of contribution from each harmonic up to the 12th for five chosen strain values (0.1, 1, 10, 40, 100% strain amplitude).The color bar indicates the extent of the contribution (purple: no contribution, brown: high contribution).

Fig. 7 .
Fig. 7. Pipkin diagram with normalized Lissajous curves from the elastic perspective for Pizza cheese cut longitudinally to the visible protein strands.Samples: L1 = Curd, L2 = intermediate sample taken during kneading and stretching, L3 = fresh cheese, L4 = cheese stored 1 week and L5 = cheese stored 2 weeks.Curves are shown for selected strains in the range from 0.1 to 150%.