Heteroscedastic sparse Gaussian process regression-based stochastic material model for plastic structural analysis

Describing the material flow stress and the associated uncertainty is essential for the plastic stochastic structural analysis. In this context, a data-driven approach-heteroscedastic sparse Gaussian process regression (HSGPR) with enhanced efficiency is introduced to model the material flow stress. Different from other machine learning approaches, e.g. artificial neural network (ANN), which only estimate the deterministic flow stress, the HSGPR model can capture the flow stress and its uncertainty simultaneously from the dataset. For validating the proposed model, the experimental data of the Al 6061 alloy is used here. Without setting a priori assumption on the mathematical expression, the proposed HSGPR-based flow stress model can produce a better prediction of the experimental stress data than the ANN model, the conventional GPR model, and Johnson Cook model at elevated temperatures. After the HSGPR-based flow stress model is implemented into finite element analysis, two numerical examples with synthetic material properties are performed to demonstrate the model’s capability in stochastic plastic structural analysis. The results have shown that with sufficient data, the distribution of the structural load carrying capacity at elevated temperatures and the variation of load–displacement curves during the loading and unloading processes can be accurately predicted by the HSGPR-based flow stress model.

The flow stress behavior takes an important role in plastic structural analysis. Many constitutive models have been proposed to describe the flow stress of material and they can mainly be divided into two types, i.e., phenomenological constitutive models and physics-based constitutive models. Phenomenological constitutive models include the Johnson Cook model 1,2 , Hockett-Sherby model 3 and Arrhenius model 4 , etc. For physical-based constitutive modelling, there are Zerilli-Armstrong model 5 and Rusinek-Klepaczko model 6 etc. The above-mentioned conventional flow stress models could achieve good description in many situations, but the applications are limited by their mathematical expression. When there is a new alloy or a complex non-linear flow behavior, the original expression will become unreliable, and a new mathematical expression often needs to be developed. For example, Samantaray et.al 7 modified the expression of the Zerilli-Armstrong model to make it suitable for the modified 9Cr-1Mo steel. However, the development of a new mathematical expression is non-trivial and can be time-consuming 8 .
To overcome the abovementioned limitations of the conventional constitutive approaches, machine learning has been introduced to model the flow stress behavior by many researchers [8][9][10][11][12][13][14] . Quan et al. 8 used one of the popular machine learning algorithms, the artificial neural network (ANN), to evaluate the flow stress of Al7050 aluminum alloy and the obtained ANN-based flow stress model showed higher accuracy than the Arrhenius model. Besides the ANN model, the support vector regression was also used by Tang et al. to predict the flow stress of Al2519 aluminum alloy 15 . The high agreement between the experimental data and the machine learningbased flow stress model was observed in the above-mentioned studies. Unlike the phenomenological constitutive models and physical-based constitutive models, machine learning-based flow stress models do not require any pre-postulated mathematical model and can predict the flow behavior more efficiently and more conveniently, as long as sufficient amount of experimental data are available 8 .
Although the machine learning based flow stress models showed good performance in the previous studies, most of the existing machine learning based models are deterministic and do not consider the material uncertainty. Due to the manufacture error of the material and the variation of the chemical compositions, even

HSGPR-based flow stress model
Basic theory of HSGPR. To achieve the sparse Gaussian process regression, the basis functions approach is used here 22,23 . As shown in the equations below, the output y is regarded as a linear combination of m nonlinear basis functions φ(x) of the inputs x: where x ∈ R 1×d is the input vector and d is the dimension of the input, y ∈ R 1 is the output, φ j (x) is the basis function. In this paper, the radial basis function kernel is chosen as the basis function for the HSGPR model and its expression is presented in Eq. (3). m is the number of the basis functions,w is the weight of basis functions and has the prior distribution N (0, A) where A is an m × m constant diagonal matrix and represents the covariance matrix of the prior distribution, p j ∈ R 1×d and Ŵ j ∈ R d×d are the parameters of the basis function and Ŵ j is a diagonal matrix. ∈ is the uncertainty of the output and is assumed to follow the normal distribution N 0, √ β 2 and β reflects the intensity of the uncertainty. To consider the input dependent uncertainty and enhancing the uncertainty estimation, the HSGPR model regards β as the linear combination of the nonlinear basis functions as given below.
where φ(x) is same as the basis functions in Eq. (2), u is the weight of the basis functions and has the prior distribution N (0, N) where N is an m × m diagonal matrix, and b is the constant part of Eq. (4) and is used to consider the situation of constant uncertainty intensity. The exponential function used here is to make β positive. In terms of training dataset X * , y * , X * = x * i n i=1 ∈ R n×d is the n × d matrix storing n training inputs vectors and y * = y * i n i=1 ∈ R n×1 is the vector storing n training outputs. The log marginal likelihood lnp y * of the training data can be expressed as www.nature.com/scientificreports/ ∈ R n×m is the n × m matrix storing the basis functions for each training inputs, B is the n × n diagonal matrix with diagonal elements B(i.i) = β x * i and represents the uncertainty intensity at the training inputs X * , and w = T B + A −1 T By * .
The parameters mentioned above can be determined by the training process. The training process is conducted by maximizing the log marginal likelihood of the training data. The optimization is achieved by using the quasi-Newton approach. After training the model, the parameters and weights in the basis functions can be updated to the optimal values. The output y at the unseen input x is normally distributed as N µ, δ 2 with where φ(x) is the updated basis functions; w , u and b are the updated parameters and = T B + A is a m × m matrix. µ and δ are the expectation and standard deviation of the prediction of the unseen input x.
HSGPR-based flow stress model. When the HSGPR model is used to model the material flow stress σ p at elevated temperatures, the plastic strain ε p and the temperature T (K) are used as the input variables and the flow stress σ p (MPa) is used as the output variable. The number of basis functions can be determined by trial and error. A larger number of basis functions can give higher accuracy but will increase the model complexity. Based on Eq. (6) to Eq. (9) and the Gaussian prior distribution assumed before, the obtained HSGPR-based stochastic flow stress model σ p ε p , T follows the Gaussian distribution with the mean of φ ε p , T w and the variance of δ 2 ε p , T , as expressed below where R is a standard normal distribution of N 0, 1 2 . The possible curve of the HSGPR based stochastic flow stress model can be sampled by sampling standard normal random variable R first and then substituting the sampled random variable into Eq. (10). The procedure of the HSGPR-based flow stress modelling is illustrated as a flow chart in Fig. 1.

Model evaluation.
To evaluate the accuracy of the obtained flow stress models, the average absolute relative error (AARE) is used here. The model with lower AARE has higher accuracy.
The expression of AARE is presented below: where σ exp pi and σ pre pi are the ith experimental flow stress and the ith predicted mean flow stress of the model respectively, and N is the number of the data.
To evaluate the uncertainty of the flow stress predicted by the models, the mean log likelihood (MLL) of the experimental data is introduced here. The MLL reflects the average likelihood of obtaining the experimental data under the given model. The model with a higher MLL means that it is more likely to reproduce the data under the given model and the material uncertainty predicted by this model has higher accuracy.
The expression of MLL is given below: where δ i is the standard deviation of the ith flow stress predicted by the model. www.nature.com/scientificreports/ Application for material modelling Material data. Without loss of generality, the experimental flow stress-plastic strain data of Al6061 aluminum alloy from the literature 28 are applied in the present work. There were six temperatures (20,100,150,200, 250, 300 °C) considered in the experiments. The experiment data obtained by the uniaxial tension test conducted on a standard "dog bone" type coupon by MTS 810 (10-kip capacity). Design of the uniaxial tension test followed the guidelines in ASTM E8. A total of 98 stress-strain curves (totally 24,011 data points) are used here. The detailed breakdown of the obtained curves is illustrated in Table 1 and all data curves are presented in Fig. 2.
Kolmogorov-Smirnov hypothesis test is conducted on the flow stress data at different temperature and different plastic strain to examine if the flow stress follows the Gaussian distribution. The p-values of the hypothesis test on the flow stress are presented in Table 2. As demonstrated in Table 2, the test does not reject the hypothesis, indicating that the flow stress obeys the Gaussian distribution, at the 5% significance level. Hence, it is reasonable to use the normal distribution for the uncertainty ∈ in the HSGPR model and Eq. (10) to describe the stochastic behavior of the flow stress.
In addition, the experimental data at temperature of 250 °C are used as the testing dataset for the HSGPR model, while the rest of data (21,523 data points) are used as the training dataset for the ANN model and HSGPR model and conventional GPR model.

Results and discussion
Four models, namely, the Johnson Cook (JC) model 29 , the ANN model, the HSGPR model, and the conventional GPR model are used to model the flow stress of the Al6061 aluminum alloys at elevated temperatures.
The Johnson-Cook flow stress model was used by many researchers to model the flow behavior of Al6061 aluminum [24][25][26][27] . The formulation of JC model is expressed in Eq. (14).  www.nature.com/scientificreports/   The parameters c 1 , c 2 , c 3 of JC model are calibrated by the stress-strain curves at room temperature (i.e. 20 °C) first. The Pearson correlation coefficients between c 1 and c 2 , and between c 1 and c 3 are -0.002 and -0.021 respectively, which means that the correlation is negligible and c 1 is independent of c 2 and c 3 30 . As the Pearson correlation coefficients between c 2 and c 3 is 0.96, c 2 and c 3 are assumed to be linearly correlated with each other here. By conducting Kolmogorov-Smirnov test, the parameters c 1 ,c 2 and c 3 follow normal distribution N 260.72, 17.32 2 MPa , N 504.26, 106.99 2 MPa and N 0.82, 0.14 2 respectively at the 5% significance level. The parameter c 4 is calibrated by the data of zero plastic strain. c 4 obeys the normal distribution N 1.277, 0.069 2 at the 5% significance level.
The HSGPR-based flow stress model for the Al6061 aluminum alloys is established by the approach described in "Section HSGPR-based flow stress model". In order to determine the suitable number of basis functions m, the relationship between the AARE on the testing dataset and m is drawn in Fig. 3. It can be observed from Fig. 3 that the improvement of the model accuracy is negligible when m is larger than 10. Hence, 10 basis functions are sufficient to produce a good estimate and more basis functions will only increase the model complexity without obvious contribution to the model accuracy. For the sake of the model efficiency and accuracy, m = 10 is adopted in the rest of the paper. With m = 10, the optimal parameters of the HSGPR model can be obtained by the training process and presented in "Appendix A" of Supplementary Information. The conventional GPR model and the ANN model are established by using the MATLAB toolbox as well for comparison. For the GPR model, the squared exponential kernel function is used as the kernel function, and the subset of the data points approximation is set as the fitting method to address the efficiency issue. The quasi-Newton method is used to optimize the parameters in the GPR model and the default value in MATLAB is used for all other settings of the GPR model. The optimal parameters of the GPR model (i.e., the signal standard deviation, the characteristic length scales and the noise standard deviation) are 228.8175, 244.4820 and 17.8537 respectively. The ANN structure is adopted from literature 9 directly. As demonstrated in the reference 9 , the network with one hidden layer consisting of 20 nodes is the suitable setting for estimating the flow stress.
Since the flow stress model needs to be called to predict the flow stress at each integration point and each increment step in the finite element analysis (FEA), the prediction time of the data-driven flow stress model will significantly influence the efficiency of the model in the FEA. The prediction time of the ANN model, HSGPR model and GPR model for 10,000 input points in the MATLAB platform are 0.019 s, 0.011 s and 26.78 s, respectively. Compared with the conventional GPR model, the HSGPR model and ANN model have lower time consumption during the prediction, which means they have higher efficiency in the FEA.
The AARE and MLL of the three machine learning models (namely, the HSGPR model, GPR model and ANN model) on the training and testing dataset are presented in Table 3. Since only deterministic prediction is given by the ANN model, only the AARE of the ANN model is calculated. The ANN model, HSGPR model and GPR model have similar AARE on the training stage, while the performance of ANN model on the testing dataset is the worst among the three models with the AARE of 12.4%. In addition, with the advanced uncertainty model, the HSGPR model presents obviously larger MLL and smaller AARE than the GPR model on the testing dataset. The larger MLL means the testing dataset is more likely to be reproduced by the HSGPR model than the GPR model.
In order to further evaluate the generalization ability of the three machine learning approaches, the crossvalidation procedure is performed here as well. The data points at one selected temperature are excluded for validation, while the rest of the data are used for training. Hence, there are totally 6 cross-validation tests. The (14)  www.nature.com/scientificreports/ AARE of the three models on the excluded validation set is demonstrated in Table 4. In the validation test 1 and 6, the validation temperature is not covered by the training dataset, so the model extrapolation ability can be tested. Notably, the HSGPR model have superior accuracy over the GPR model and ANN model on extrapolation. Besides, the performance of the HSGPR model on the interpolation validation set (No. [2][3][4][5] is also better than the other two models. It can be concluded that the HSGPR model has better generalization ability than the conventional GPR model and ANN model on flow stress prediction. As a conventional constitutive model, the JC model is directly evaluated using the whole dataset (including all the experimental data). For comparison, the three machine learning approaches are also evaluated on the whole dataset and the results are listed on Table 5. Due to the limitation of the mathematical expression, the performance of the JC model is worse than the machine learning approach. The HSGPR model has the smallest AARE and the largest MLL among all the models, which means that the HSGPR model shows better performance than the other models.
As for the material uncertainty predicted by the models, it can be visualized by the 95% confidence region (i.e. mean ± 2 × standard deviation ) and the coefficient of variation (COV, i.e. standard deviation/mean ) of the flow stress model and is plotted in Figs. 4, 5, 6, 7. Since the ANN cannot provide the material uncertainty, only the HSGPR-based flow stress model, GPR-based flow stress model and JC model are plotted and discussed in the next.
At the low temperatures (20 °C to 150 °C), it can be observed from Figs. 4, 5, 6, 7 that the 95% confidence region given by the three models could cover most of the data at low temperatures. Besides, as shown in Fig. 4, the COVs given by the three models are similar and are close to the COV given by the experimental data at the temperatures of 20 °C to 150 °C.
At the high temperatures (200 °C to 300 °C), the mean curves and material uncertainty estimated by the HSGPR model are consistent with the experimental data, while the material uncertainty is overestimated by the GPR model and the COV given by GPR model is much larger than the experimental value. It is notable to mention that the advanced uncertainty model within HSGPR model results in the better uncertainty estimation  www.nature.com/scientificreports/ than the GPR model. In addition, as limited by its mathematic expression, the JC model cannot describe the severe temperature softening of Al6061 alloys at high temperatures. As shown in Fig. 7, the estimated flow stress given by the JC model deviates from the experiment data at high temperature, which causes the high AARE of the JC model in Table 5. The poor performance of the JC model for the Al6061 alloy at high temperatures is consistent with the results in the work of Fan et.al. 25 . Although it is possible to modify the JC model to improve the accuracy, those modifications are material-dependent and need to postulate the mathematical expression as well 31 . Establishing a new mathematic expression of the modified JC model for each new material is always a time-consuming work and sometimes a suitable modification is hard to find. Different from the JC model, the HSGPR-based flow stress model is a material-independent data-driven flow stress model. Without postulating the specific mathematical form, the HSGPR framework can be used to establish a flow stress model directly from the material data with high accuracy. In addition, it should be admitted that the current HSGPR based model use the uncertainty model ∈ with the normal distribution to approximate the uncertainty of material data. If the data deviated from the normal distribution, the uncertainty ∈ with other distribution needs to be used. The new distribution of ∈ can be selected by trial and error, and the distribution, which can produce the maximum MLL, is regarded as the suitable one.  www.nature.com/scientificreports/ where σ is the equivalent von Mises stress, σ p is the HSGPR-based flow stress, ε p is the von Mises equivalent strain, T is temperature, and [σ 11 , σ 22 , σ 33 , σ 12 , σ 13 .σ 23 ] is the Cauchy stress vector.
(15) f σ , ε p , T = σ − σ p ε p , T www.nature.com/scientificreports/ The yield function with the HSGPR-based flow stress model is called at each increment of the FEA to check the yield state of each element integration point and update the flow stresses of the element integration points which have reached plasticity. In order to update the flow stresses, the stiffness of the HSGPR-based flow stress model, that is the gradient of the flow stress model σ p with respect to the equivalent plastic strain ε p and temperature T , respectively, is queried in the FEA 32 . The corresponding stiffness (i.e. ∂σ p /∂ε p and ∂σ p /∂T ) of the HSGPR-based flow stress model is expressed as: Stochastic plastic structural analysis procedure. The Monte Carlo Simulation (MCS) method is applied to achieve the stochastic plastic structural analysis using the HSGPR-based model. The possible flow stress relations estimated by the HSGPR model are generated by following the flow chart in Fig. 1 and samples of other possible basic random variables, such as elastic modulus (E), can also be generated according to their probability distribution. For each sample, a finite element model is set up and a complete analysis is performed using Abaqus. The procedure of performing the plastic stochastic structural analysis with the HSGPR-based model is presented schematically in Fig. 8 (20,100,150,200,250, 300 °C) are considered in the training dataset. For each temperature, 20 stress-strain curves (24,000 data points), 100 stress-strain curves (120,000 data points), and 200 stress-strain curves (240,000 data points) are generated for studying the effect of the dataset size. Hence, there are totally nine synthetic flow stress datasets with three uncertainty levels and three dataset sizes generated for training the HSGPR model. The Poisson ratio of the material is set as 0.33 and the elastic modulus is set as a temperature-dependent normal random variable 35 and is expressed as www.nature.com/scientificreports/ where T (K) is the temperature, T r is the room temperature 293.15 K and the room temperature elastic modulus E 0 is a normal random variable with a mean of 69,500 MPa and a standard deviation of 8,548.5 MPa 35 . After the training process, the HSGPR-based flow stress model is implemented into the FEA and the model is validated by the beam problem adopted from 36 . The geometry, boundary conditions and loading conditions of the beam are depicted in Fig. 9 and the thickness of the beam is 35 mm. The failure strain is set as 0.1 16 . Except for the elastic modulus and flow stress, other uncertainty sources are not taken into account in the current validation process. Abaqus is used to conduct the finite element analysis. The displacement-controlled loading is applied at the middle point of the beam model. The mesh convergence study is conducted with all the parameters being set to their mean values. The linear 4-node doubly curved shell element (S4R) is used in the Abaqus model. As depicted in Fig. 10, four different mesh size strategies are studied in the present work. The load-displacement curve of each case is presented in Fig. 11. The relative errors of the load capacity (i.e. peak load) of the cases with 10 mm, 5 mm and 1 mm mesh relative to the case with 0.5 mm mesh are 1.62%, 0.35% and 0.32% respectively. The case with 5 mm mesh, which could provide sufficient accuracy for estimating the load capacity and acceptable computational efficiency, is used in the current study. The mesh distribution of the beam is presented in Fig. 9.
In order to determine the number of MCS, the stochastic structural capacity of the beam is evaluated by the HSGPR model at the different number of iterations. The expectation and standard deviation of the structural capacity with respect to the number of MCS is presented in Fig. 12. It can be noted from the figure that both the expectation and standard deviation given by the MCS converge at the number of 1000. Hence, the number of MCS is set as 1000 for the structure analysis.
The probability density distributions of the structural load capacities given by the HSGPR-based flow stress model and the targe model at the room temperature and the elevated temperatures are presented in Fig. 13. The corresponding mean and standard deviation of the load capacities are illustrated in Tables 6 and 7, respectively. It can be observed from Table 6 that the expected capacity given by the HSGPR-based flow stress model have high accuracy with the relative error smaller than 4% even under the high uncertainty level and small dataset size. When the dataset size is large enough, the influence of increasing dataset has little effect on the accuracy of expected capacity. The slight error variation may be due to the model itself. For the variability of the load capacity, the probability density distribution obtained by the HSGPR model coincides with the distribution obtained by the target model at the low uncertainty level (COV = 0.05), as shown in Fig. 13. For the high uncertainty level, the HSGPR model requires larger dataset size to achieve good estimation of the COV of load capacity. With the sufficient material data, the relative error of the standard deviation given by the HSGPR model is smaller  www.nature.com/scientificreports/ than 3%. As presented in Tables 6 and 7, the effect of temperature on the accuracy of the HSGPR model is not obvious at all uncertainty levels, which means that the HSGPR-based flow stress model is suitable for the stochastic plastic structural analysis at elevated temperatures with sufficient dataset. In addition, it is worth noting that the HSGPR model can give good estimation on the load capacity at the temperature of 175 °C which is not considered in the synthetic dataset.
Punch problem. To evaluate the performance of the HSGPR-based model in the loading and unloading process, the punch problem given by Huang et. al. 37 is used as the benchmark test. The exponential flow stress law σ p = e 1 + e 1 (e 2 + ε p ) e 3 is set as the target model. To consider the material uncertainty, the parameters e 1 , e 2 and e 3 are set as the normal random variables with the mean of 0.05 MPa, 0.00002 and 0.3 respectively and the COV of 0.05. A total of 200 stress-strain curves (40,000 data points) at the room temperature are generated for training the HSGPR model. The Poison ratio is set as 0.33 and the elastic modulus is set as the normal random variable with the mean of 1 MPa and the standard deviation of 0.05 MPa. The HSGPR-based flow stress model obtained after the training process is evaluated by the punch problem 37 . The geometry of the problem is depicted in Fig. 14. The left side and the top right corner of the block are only fix in the horizontal direction and the bottom of the block is only fix in the vertical direction. During the loading stage, the vertical displacement of the top of the block u 0 is gradually increased to 0.07 mm, which can be regarded as a quasi-static process. During the unloading stage, the direction of the displacement is reversed until u 0 = 0.02 mm. Similar to Huang et. al. 37 , the block is divided into 100 quadratic quadrilateral elements for the FEA analysis. The load-displacement curve of the node A is recorded.  www.nature.com/scientificreports/ For the deterministic analysis, the obtained load-displacement curve is plotted in Fig. 15. The load displacement simulated by using the feed forward neural network (FNN)-based plasticity model by Huang et al. 37 is also presented in Fig. 15. It can be observed from Fig. 15 that the HSGPR model and FNN model given by Huang et al. 37 have the similar performance with both being close to the target model during the loading and unloading processes. And the average relative error of the HSGPR model is 0.75%. However, unlike the HSGPR model, the FNN model cannot capture the material uncertainty and be used in the stochastic plastic structural analysis.
Similar to the beam problem, 1000 Monte Carlo simulations are sufficient to provide the converged estimation for the punch problem as well. By conducting 1000 Monte Carlo simulations using the HSGPR model, the 95% confidence region (i.e. mean ± 2 × standard deviation ) of the load-displacement curve can be obtained and is plotted in Fig. 16. As depicted in Fig. 16, the confidence region given by the HSGPR model is in agreement with the target model and the variation of the load-displacement curve is depicted accurately by the HSGPR model in both the loading and unloading stages.
In order to investigate the influence of material uncertainty on the plastic deformation, the equivalent plastic strain of the block after loading is analyzed by using the HSGPR model. Two situations corresponding to the upper bound and the lower bound, as labeled in Fig. 16, are shown in Fig. 17. It can be observed from Fig. 17 that the plastic deformation of the upper bound situation in the lower left conner of the block is smaller than www.nature.com/scientificreports/  www.nature.com/scientificreports/ the lower bound situation. However, both the upper bound and the lower bound are the possible situations of the structural responses. The variation of the deformation is caused by the uncertainty inherent in the material. Hence, it can be concluded that the deformation of the structure cannot be described completely without sufficiently considering the flow stress uncertainty.

Conclusions
Understanding the uncertainty of material flow behavior is important in the stochastic plastic structural analysis. In this paper, a data-driven approach, the heteroscedastic sparse Gaussian process regression (HSGPR)-based flow stress model, is proposed to capture the flow stress behavior and the associated uncertainty directly from the available material data. The established HSGPR-based flow stress model is verified by the experimental data of the Al6061 aluminum alloy and is compared with the conventional GPR model, the ANN model and the Johnson Cook (JC) model. The flow stress of the Al 6061 aluminum alloy at elevated temperatures predicted by the HSGPR model are more accurate than those by the ANN model, the GPR model and the JC model, with the AARE of 6.5% and the MLL of -3.74 on the testing dataset. Besides, unlike the ANN model, which provides only the deterministic flow stress, the flow stress uncertainty is successfully estimated by the HSGPR based model as well.
The HSGPR-based stochastic flow stress model is then implemented into finite element method for stochastic plastic structural analysis. Two numerical examples with synthetic flow stress data, that is the beam problem and the punch problem, are used to verify the accuracy of the HSGPR model in stochastic plastic structural analysis. With the sufficient material data, the expectation and standard deviation of the structural capacity at elevated temperatures have the relative errors less than 1% and 3% respectively in the beam problem. In the punch problem, the variation of the loading and unloading paths given by the HSGPR model is consistent with that given  www.nature.com/scientificreports/ by the target model. The influence of material uncertainty on the structure plastic deformation is successfully identified by the HSGPR model.

Data availability
The raw/processed data required to reproduce these findings cannot be shared at this time due to technical or time limitations. However, these data will be shared upon request to the corresponding author.