Isaac Y.
Miranda-Valdez
*,
Aaro
Niinistö
,
Tero
Mäkinen
,
Juha
Lejon
,
Juha
Koivisto
and
Mikko J.
Alava
Department of Applied Physics, Aalto University, P. O. Box 15600, 00076 Aalto, Espoo, Finland. E-mail: isaac.mirandavaldez@aalto.fi
First published on 20th March 2025
Mathematical modeling is a powerful tool in rheology, and we present pyRheo, an open-source package for Python designed to streamline the analysis of creep, stress relaxation, small amplitude oscillatory shear, and steady shear flow tests. pyRheo contains a comprehensive selection of viscoelastic models, including fractional order approaches. It integrates model selection and fitting features and employs machine intelligence to suggest a model to describe a given dataset. The package fits the suggested model or one chosen by the user. An advantage of using pyRheo is that it addresses challenges associated with sensitivity to initial guesses in parameter optimization. It allows the user to iteratively search for the best initial guesses, avoiding convergence to local minima. We discuss the capabilities of pyRheo and compare them to other tools for rheological modeling of soft matter. We demonstrate that pyRheo significantly reduces the computation time required to fit high-performance viscoelastic models.
Mathematical modeling in rheology is a gateway to understanding the structure–property relationship for soft matter. However, choosing a model, data processing of the experimental measurements, and curve-fitting routines can represent a steep learning curve in conducting and interpreting viscoelastic experiments, due to the highly nonlinear nature of the behavior. For example, rheological models based on fractional order derivatives commonly require a curve-fitting routine that involves computationally expensive operations, such as the Mittag–Leffler function—most commonly represented by an infinite sum of terms containing the gamma function.8–10 Another challenge in the mathematical modeling of viscoelasticity is defining a cost function that allows selecting a model and further finding its parameters.11–13 Therefore, choosing a model and inferring its parameters are two critical decisions that are highly uncertain and have non-unique solutions.13
Currently, there are numerous open-source rheological tools. For example, Boudara et al.14 developed RepTate which offers comprehensive tools for analyzing linear and nonlinear rheological data. In particular, RepTate provides an interface to analyze entangled polymer melts using theoretical frameworks such as the Likhtman–McLeish theory and multi-mode Maxwell analysis. Other remarkable efforts include the work of Luciano et al.,15 who designed oreo to analyze nonlinear rheological data, and the work of Tassieri et al.16 who developed i-Rheo to infer linear viscoelastic properties with Fourier transform analysis. Additionally, Shanbhag17 developed pyReSpect as a tool to extract relaxation spectra from stress relaxation experiments. Nonetheless, despite these great efforts, no open-source tools for Python integrate rheological frameworks based on fractional calculus. Therefore, this work introduces pyRheo, an open-source Python package that assists in model selection and fitting procedures for several types of rheology experiments.
pyRheo focuses on streamlining the mathematical modeling of rheological data obtained in the linear viscoelastic regime using fractional rheology.18–20 First, pyRheo uses machine intelligence to suggest a rheological model likely to describe a provided dataset. Then, it allows users to fit the proposed model or choose a different one. pyRheo enables the user to automatically or manually choose from several fractional rheological models. pyRheo focuses on fractional order viscoelastic models which have proved to be able to provide valuable insights into soft materials, such as food gels, structured fluids, biological tissues, cells, and polymers.21–27 For example, fractional models have allowed linking the microstructure of colloidal gels to their relaxation spectrum and elasticity.19,20
pyRheo is highlighted as a tool for analyzing rheological data readily and as an interface that can be coupled with machine learning algorithms widely available for Python.28 The following sections demonstrate pyRheo's features. We present a set of robust validations against experimental results and other public toolkits to check the accuracy and computational performance of our code package in characterizing soft materials such as biological tissues, polymers, foams, foodstuff, and gels. We note that pyRheo is available via its GitHub repository, where all the Python scripts to compute every single example presented in this paper and its ESI† are available as Jupyter Notebooks. Furthermore, we have created a simple graphical user interface (GUI) for those users whose programming skills may limit their access to pyRheo. The GUI file can be found in the GitHub repository.
First, we test the performance of pyRheo using the creep data measured for a perihepatic abscess reported by Shih et al.30 and the stress relaxation data of a fish muscle reported by Song et al.2. The first step in pyRheo's workflow is to import the rheological data of creep compliance J(t) and relaxation modulus G(t) into the MLP model. In the cases shown in Fig. 2a and b, the MLP model classifies the data from creep as FractionalKelvinVoigt and stress relaxation as FractionalMaxwell. FractionalKelvinVoigt consists of two springpots connected in parallel, whereas FractionalMaxwell is built by two springpots connected in series (see insets in Fig. 2).2 The predicted model is automatically fitted to each dataset by pyRheo, and the fitting results are depicted by the dashed lines in Fig. 2a and b.
![]() | ||
Fig. 2 Creep and stress relaxation data classified and fitted with pyRheo. (a) Creep compliance J(t) of a perihepatic abscess sample. The curve is fitted using the auto function in pyRheo, which classifies the data as a FractionalKelvinVoigt. (b) Relaxation modulus G(t) of a fish muscle classified and fitted with FractionalMaxwell. The raw data of the perihepatic abscess was reproduced from Shih et al.30 and the data of the fish muscle from Song et al.2 |
In Fig. 3, we showcase an instance of coupling Lennon et al.'s29 package with pyRheo to analyze the small amplitude oscillatory shear (SAOS) data of an interpenetrating-network hydrogel made of cellulose nanofibers and methylcellulose. We import the master curve data from G′(ω) and G′′(ω) to the MLP classifier of pyRheo. The MLP classifies the G′(ω) and G′′(ω) data as a FractionalKelvinVoigt. The result of fitting this model to the master curve data is depicted by the solid and dashed lines in Fig. 3a. Furthermore, we demonstrate in Fig. 3b how to utilize the model predictions generated with pyRheo to easily construct other visualizations such as Cole–Cole diagrams. This is feasible because pyRheo stores the model results as an object the user can call to, for example, predict the material response according to a specified ω range. Consequently, this flexibility allows for estimating model predictions to higher and lower ω values.
pyRheo's model fitting and prediction tools are specifically designed to address fractional order viscoelastic models, as these models provide a more succinct representation of viscoelastic phenomena compared to traditional multi-mode fitting methodologies, such as generalized Maxwell models. This advantage becomes particularly pronounced when characterizing the linear viscoelastic response of entangled polymer melts.
As illustrated in Fig. 4, we compare the fitting outcomes for a polyisoprene melt using the FractionalZener model (depicted in Fig. 4a) with those obtained from a generalized Maxwell model employing eight modes (shown in Fig. 4b). Both models effectively capture the behavior of the storage modulus G′(ω) and the loss modulus G′′(ω). Notably, the fractional order model necessitates only six parameters, whereas the generalized Maxwell model requires 16 parameters, thereby illustrating the enhanced efficiency of the fractional framework. Schmidt et al.31 have thoroughly discussed the potential of fractional frameworks in comparison to generalized approaches.
![]() | ||
Fig. 4 Comparison between pyRheo and RepTate frameworks applied to describe the small amplitude oscillatory shear (SAOS) data of a polyisoprene melt (data from Boudara et al.14). (a) pyRheo fitting using FractionalZener. (b) RepTate fitting using a generalized Maxwell model with eight modes. |
This section is complemented by additional results in ESI Note 4,† where we present more fitting routines, which include materials such as mucus, foams, polymer networks, gels, plastics, food colloids, and polysaccharides. The examples included in ESI Note 4† present data analysis from steady shear flow experiments, which are not detailed in the main article due to their lower computational complexity. For transparency, every demonstration with pyRheo is available as a Jupyter Notebook on the pyRheo GitHub page, which users can test and adapt to suit their needs.
![]() | (1) |
There are several methods for numerically computing the ML function, either through the numerical inversion of its Laplace transform or by using mixed techniques, including Taylor series, asymptotic series, and integral representations.32 Notable examples of these methods can be found in the algorithms developed by Garrappa33 and Podlubny.34
In fitting routines, the computational demands of the ML function can become increasingly sensitive to the size of the dataset. This often leads to exponential growth in computation time as the dataset expands. The latter is common in master curve fitting and creep and stress relaxation tests, where the sampling rate is higher than in SAOS and steady shear flow tests. Current methods for reducing the computation time of the ML function typically involve downsampling. However, this process introduces uncertainty and can lead to non-unique outcomes in the fitting. The computational demand of the ML function is also evident when the fitting involves iterative optimization of multiple parameters; poor initial parameter guesses can result in slower convergence or even lead to convergence to a local minima. Consequently, finding the best model parameters may require restarting the optimization with a different initial guess for the model parameters.
To the best of our knowledge, using Padé approximations to compute rheological models has not been extensively researched. Thus, pyRheo exploits the Padé approach to reduce the computation time spent in fitting fractional rheological models by implementing the ML function based on the global Padé approximations proposed by Zeng and Chen32 and Sarumi et al.35. In Fig. 5, we compare the performance of pyRheo against popular methodologies for fitting fractional rheological models, which are based on the MATLAB toolkit published by Song et al.2 and the RHEOS package for Julia programming language.36 The MATLAB toolkit uses Garrappa's algorithm33 to compute the ML function, whereas RHEOS uses Gorenflo et al.'s37 approach.
![]() | ||
Fig. 5 Performance comparison of fitting routines across pyRheo (Python), MATLAB, and RHEOS (Julia) for the relaxation modulus G(t) of polyethylene (PE) and fish muscle. (a) Fitting results for PE from each implementation using FractionalZenerSolidS. (b) Normalized computation times tn for fitting routines on PE; RHEOS required downsampling to 20% of the dataset, while pyRheo and MATLAB processed the full dataset. (c) Fitting results for fish muscle from each implementation using FractionalMaxwell, based on Song et al.2. (d) Normalized computation times tn for fish muscle fitting; RHEOS required downsampling to 10%, while others used the complete dataset. All computations were on a system with an Intel Core i5-12600K CPU at 3.7 GHz, 31 GB RAM, and 1 TB SSD, running Ubuntu 20.04. |
Accordingly, in Fig. 5a, we fitted a FractionalZenerSolidS to the relaxation modulus G(t) of a polyethylene (PE) sample to show the applications of pyRheo in other fields of soft matter. For the examples reported here, we chose Pade32 in pyRheo; in other words, a second-order global Padé approximation. Fig. 5b displays the computation time spent in fitting the stress relaxation data of PE. The fitting requires computing the one-parameter ML function (i.e., b = 1). pyRheo, MATLAB, and RHEOS yield similar parameter values, as seen in Table 1. However, it is essential to note that the computation times tn vary significantly among the three implementations. For the stress relaxation data of PE, pyRheo identifies the optimal parameters one to three orders of magnitude faster than MATLAB and RHEOS, respectively. We normalize the computation times by the size of the dataset to enable a fair performance comparison across the different implementations. In the case of RHEOS, we downsampled the PE dataset to contain 20% (tn is obtained by multiplying the real time by a factor of five) of the original one since the computation time extended beyond the capabilities of a desktop computer.
In Fig. 5c, we show again the stress relaxation data of the fish muscle from Fig. 2b. Fig. 5d shows that the computation of FractionalMaxwell is more time-consuming than that of the FractionalZenerSolidS. The latter is due to the multi-parametric nature of the ML function used in the FractionalMaxwell. Again, in Table 1, we observe that the three implementations find similar parameter values for the fish muscle. Nonetheless, the computation times with pyRheo are shorter than those needed by MATLAB and RHEOS. pyRheo leverages the computational efficiency of the ML function thanks to the global Padé approximation. Again, we normalized the computation times for the fish muscle data to enable a fair performance comparison across the different implementations. In the case of RHEOS, we downsampled the datasets to contain 10% of the original ones (i.e., tn is the real computation time scaled by a factor of ten).
In addition to the global Padé approximation, we also programmed in pyRheo the option to use Garrappa's algorithm33 for the evaluation of the Mittag–Leffler function. The algorithm was adapted from its MATLAB script.38 The inclusion of Garrappa's algorithm allows users to benefit from its robust computation method, especially in cases where the Mittag–Leffler function needs to be evaluated for parameters that pose challenges for the global Padé approximation. This flexibility is crucial as it ensures that accurate and reliable results can be obtained across a broader range of applications and parameter settings, reinforcing pyRheo's position as a trusted tool for researchers and engineers in rheology.
In ESI Notes 5 and 6,† we detail the global Padé approximations and provide examples demonstrating the differences between the global Padé approximation and Garrappa's algorithm when computing the FractionalMaxwellGel, FractionalMaxwellLiquid, and FractionalMaxwell. These examples showcase the accuracy and reliability of each algorithm in various scenarios. Such detailed comparisons help users make informed decisions about which algorithm to employ for their specific needs, thereby enhancing the overall utility and effectiveness of pyRheo in addressing complex rheological analyses.
Our work demonstrates how to reduce the computational cost of fitting routines involving the ML function, which is part of many models constitutive equations for creep compliance J(t) and relaxation modulus G(t). We decreased the computation times by utilizing Padé approximations to compute the ML function. This advancement allows for exploring a wider range of models and methodologies, requiring fewer resources and less time. Besides faster computation of models with the ML function, the Padé approximation enables pyRheo to offer a solution to a common problem in fitting routines of rheological models, which is sensitivity to the initial guesses. Commonly, a bad choice of initial guess might lead the parameter optimization process to converge to local minima. pyRheo presents two solutions to determine the initial guesses: random search and Bayesian optimization (BO).28
In the data analyzed in Fig. 3 and 5, we have implemented a random search of initial guesses. In all cases, we have fixed a maximum of ten restarts of the optimization algorithm that seeks to minimize the weighted residual sum of squares RSSwi. As shown in Fig. 5, our brute-force approach combined with the global optimization algorithms from SciPy39 yields faster computations than the MATLAB and RHEOS implementations. On the other hand, in ESI Note 7,† our BO approach was also shown to be effective in finding suitable initial guesses for fitting the creep data of a mucus gel. Our work reveals how random search and BO methodologies, techniques that are used in hyperparameter tuning of machine learning, can be adapted to traditional fitting routines.
The fitting tools available in pyRheo, along with its integrated machine intelligence, position it as a valuable resource for developing automated laboratories capable of conducting high-throughput testing and analysis. In the future, the Multi-Layer Perceptron (MLP) classifier could assist high-throughput rheometers in reformulating and optimizing materials. However, it is important to note that the current MLP classifier integrated into pyRheo cannot label rheological data related to Zener models. This limitation arises because the responses of Zener models overlap with those of Maxwell and Kelvin–Voigt models. One potential strategy to address this issue would be to use multi-label classifiers.
We have also integrated in pyRheo an option for the user to simply call a model's function (evaluators class). This option allows the user to compute the fractional models by assigning fixed values to the model parameters, as shown in the plots in ESI Note 2.† As shown by Miranda-Valdez et al.,40 this evaluator class offers the user flexibility to design their own problems.
The accuracy (with synthetic data) of the MLP classifiers ranges from 70 to 80%. We suggest using the auto method as a first approximation to identify the type of rheological behavior. More detailed information about the machine learning training process and performance is disclosed in the ESI.†
![]() | (2) |
Users may define their own initial guesses and parameter bounds (automatic bounds and random initial guesses by default). Then, to minimize RSSwi, users can choose from several minimization algorithms implemented on SciPy,39 such as Nelder–Mead, Powell, and L-BFGS-B (Powell by default).
As shown in Table 2, an advantage of using pyRheo is that it addresses the challenges associated with sensitivity to initial guesses in parameter optimization. In other words, if an initial guess is close to a local minimum, the minimization algorithm may converge there instead of the global minimum. Therefore, pyRheo allows the user to restart the fitting process multiple times with random initial parameter values and then take as the final result the iteration with the lowest RSSwi. By generating a diverse set of random starting points, this brute-force approach increases the likelihood of exploring different regions of the parameter space, thus avoiding local minima.
Another method that pyRheo offers for defining initial guesses is Bayesian Optimization (BO).41–43 In this approach, pyRheo creates a mapping from the parameter space to the error space
using Gaussian Process Regression (GPR), represented as
where g is the Gaussian Process. The surrogate model
is developed by computing the constitutive equation of the target model with fixed parameter values and then recording the difference (residuals) between this computation and the data being analyzed. The goal of BO is to minimize ε by exploring various combinations of parameter values, guided by an acquisition function known as expected improvement, which balances exploration and exploitation of the parameter space. Afterward, pyRheo uses the BO solution as the initial guess for the minimization algorithm. Miranda-Valdez et al.40 present this methodology in detail and expand it further.
Footnote |
† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d5dd00021a |
This journal is © The Royal Society of Chemistry 2025 |