Modeling resilient modulus of subgrade soils using LSSVM optimized with swarm intelligence algorithms

Resilient modulus (Mr) of subgrade soils is one of the crucial inputs in pavement structural design methods. However, the spatial variability of soil properties and the nature of test protocols, the laboratory determination of Mr has become inexpedient. This paper aims to design an accurate soft computing technique for the prediction of Mr of subgrade soils using the hybrid least square support vector machine (LSSVM) approaches. Six swarm intelligence algorithms, namely particle swarm optimization (PSO), grey wolf optimizer (GWO), symbiotic organisms search (SOS), salp swarm algorithm (SSA), slime mould algorithm (SMA), and Harris hawks optimization (HHO) have been applied and compared to optimize the LSSVM parameters. For this purpose, a literature dataset (891 datasets) of different types of soils has been used to design and evaluate the proposed models. The input variables in all of the proposed models included confining stress, deviator stress, unconfined compressive strength, degree of soil saturation, soil moisture content, optimum moisture content, plasticity index, liquid limit, and percent of soil particles (P #200). The accuracy of the proposed models was assessed by comparing the predicted with the observed of Mr values with respect to different statistical analyses, i.e., root means square error (RMSE) and determination coefficient (R2). For modeling the Mr of subgrade soils, percent passing No. 200 sieve, optimum moisture content, and unconfined compressive strength were found to be the most significant variables. It is observed that the performance of LSSVM-GWO, LSSVM-SOS, and LSSVM-SSA outperforms other models in predicting accurate values of Mr. The (RMSE and R2) of the LSSVM-GWO, LSSVM-SSA, and LSSVM-SOS are (6.79 MPa and 0.940), (6.78 MPa and 0.940), and (6.72 MPa and 0.942), respectively, and hence, LSSVM-SOS can be used for high estimating accuracy of Mr of subgrade soils.

Structural responses of pavement are crucial requirements to assess the quality of pavement construction materials under traffic loading changes 1,2 . Thus, based on the Association of State Highway Officials (AASHTO), the pavement characteristics should be qualified in the road design and used 2,3 . Resilient modulus (Mr) of pavement materials is one of the important characteristics, which is defined as the ratio of dynamic deviatoric stress to the recoverable strain under acyclic pulse load 2,4 . For subgrade soil and unbound granular materials, Mr is used to measure the elastic modulus of soil layers at a given stress level, and describe the non-linear stress-strain of soils under dynamic loads 2 . Many studies were devoted to investigating the effect of different factors on the Mr of subgrade soils. These studies concluded that stress state, dry density, aggregate gradation, amount of fines (materials passing the standard US sieve No. 200), moisture content/matric suction, particle shape, and aggregate type have a significant impact on the Mr 3,5-9 .
Normally, Mr determination uses expensive experimental testing and is time-consuming; furthermore, the spatial variability of soil properties and the nature of the test protocols have made Mr's determination complex and inexpedient. Therefore, many researchers utilized the machine learning technique to estimate the Mr of pavement materials 3,4,[10][11][12][13] . Predicting the Mr of subgrade soils is important in building a relationship between input and output variables for modeling those values. Consequently, the accurate modeling of Mr can help in

Least square support vector machine (LSSVM)
In order to improve the performance of SVM, the LSSVM was proposed by Suykens and Vandewalle 40 . LSSVM methods work out linear matrix problems with fewer constraint conditions 41,42 . The main advantages of LSSVM are that it overcomes the SVM drawbacks in computation cost and uncertainties in structural parameter determination 41 . Compared to SVM, LSSVM is a more powerful computation in solving nonlinear and smalldata problems 43 . LSSVM is used for classification and regression problems. This study aims to develop a highperformance hybrid machine learning model using a regression-based machine learning model, i.e., LSSVM. The regression modeling of LSSVM can be summarized as follows 35,36 : For training of given points l, x i , y i |i = 1, 2, 3, . . . l , x i ∈ R n is input variables and y ∈ R is the output variable; the regression fitting output of LSSVM can be expressed as: Based-optimizing formulas: where, w, b and β are weight vector, deviation, and error variation, respectively; ϕ(.) denotes the mapping function, and c ∈ R + is penalty parameter. Largrange method is utilized to solve the above equations 40 .
In this study, the radial-basis function is selected, which was performed in similar studies and found to be the best in modeling non-linear behaviors 11,35 , and it can be expressed as follows 36 : The σ is the kernel function width. Therefore, the fitting model of the final output can be expressed as: Here, the main disadvantage of LSSVM is that accuracy depends on the regularization parameter ( γ ) and kernel function parameter, i.e., kernel width ( σ ) to improve the modeling results. Although the effectiveness of reconstructed input datasets with optimal parameters of conventional LSSVM in some cases 43 , it may have some inherent bias with changing the trend of other cases. In addition, time consuming and a priori knowledge requirements may decrease the model accuracy. For that swarm, algorithms have been developed to provide effective LSSVM parameters that reduce the bias of changing the data inputs to overcome the time consuming and a priori information requirements. This study uses the following swarm intelligence meta-heuristics algorithms to optimize those parameters. Herein, it should be mentioned that the proposed model can be used to optimize Mr within the limitation of input datasets.

Swarm intelligence meta-heuristic algorithms
Particle swarm optimization (PSO). PSO is a population solution-based social search behaviour of swarm members "particles" 35 . It begins with random initialization of particles in the search space to build their own and neighbours' previous successful attempts 35,38 . This aims to find the best position of particles through change their location and updating their velocity in the research space 20,38 . Mathematically, the particles X i = (x i1 , x i2 , . . . ., x iD ) and these situations with best fitting can be represented in the best current and global positions, these are P i = (p i1 , p i2 , . . . , p iD ) and P g = (p g1 , p g2 , . . . , p gD ) . These can be attained through the best fitting function Pbest , and the best global function Gbest , respectively 38 . The velocity of particles is represented The following equations represent the velocity and position of particles updating in each iteration 38 : where, r 1 and r 2 are random values ǫ(0, 1) , c 1 and c 2 are acceleration coefficients, and ω is the inertia-weight factor. Grey wolf optimization (GWO). GWO is also a population solution that simulates the social behaviour of a grey wolf pack 44 . In GWO, four categories are divided, are alpha, beta, delta and omega 36 . The leaders of the whole pack represent as alpha category (higher category). Beta group helps the leader's group in implementing commands on other lower categories. Delta group uses to fulfil above commands and controls omega . The omega group mainly follows all leaders' commands by superior departments. The hunting plan contains three steps, that are identifying and chasing the prey, encircling and harassing prey until it stops resilience, and attacking on prey 36,44 .
The GW population is assumed in the optimization process as n with unknown d-dimensional search space 36 . GW positions can be expressed as The best fitting solution of alpha, beta, and delta can be considered as X α , X β , and X δ , respectively.
The hunting process, mathematically, can be modelled as follows: where, X p and X w are the position of prey and wolf, respectively; and A and C are coefficients and can be calculated as: where, � r 1 and � r 2 ∈ (0, 1) ; a is a coefficient decreases linearly from 0 to 2 with increasing the iteration numbers.
The best positions for the best solution in this optimization method are determined based on the hunting process. The alpha guidance performs in the hunting process. Beta and delta might follow alpha in trapping a prey to find the best solution for the prey. Then the omega are pressed to follow and update positions according to the best positions of alpha, beta, and delta . Mathematically, the positions can be expressed as follows 36 :  45 . The three symbiotic interactions can be summarized as follows 23,45 : (a) The interactions of an organism with another organism are commonly beneficial in the symbiotic mutualism stage. This stage can be expressed as follows: where, x i andx ii represent the ith and iith organism vectors of the ecosystems in i and ii, i = ii . x best denotes the best organism in the current iteration, x inew and x iinew denote the respective organism for x i andx ii after their interaction, and f is the fitness function.
(b) In the commensalism stage, the interactions of an organism with another one benefit that organism and possess no effect on the other organism. This can be mathematically expressed as follows: (c) In the parasitism phase, the interactions of an organism with another one benefit that organism and harm the other organism, and this can be expressed as: where, x parasite is the artificial parasite organism generated to compete with x ii ; LB and UB represent the lower and upper bound of the problem.
Salp swarm algorithm (SSA). SSA is a new population solution that simulates the behaviour of slaps in oceans throughout locomotion for finding the best solution to optimization problems 46,47 . The details of this method can be found in 48 . The summary of this method can be proposed as follows 37 : slaps divided into leaders and followers to mimic the best position solution. The followers follow the leader in the chain. The best solution can be found through the leader who is the front in the chain using the model up to a particular iteration. For n problem variables, the position of the slaps is saved in the n-dimension search space. This can be stored in a two-dimension matrix k , which presents the position of the leader and it can be expressed as follows: where, P represents the food source, and P j is the position of food in the jth dimension. u j andl j denote the upper and lower bounds of jth dimensions, respectively. d 2 andd 3 ∈ (0, 1) , and d 1 represents a coefficient can be determined as: where t and T are the current and total number iterations, respectively. Here, the updating of the follower's position after each iteration can be calculated based on Newton's law of motion as: where k i j represents the ith follower particle in the jth dimension; where, i = 1 for the leader position (Eq. 16), and i ≥ 2 is for the followers. t is time, v 0 denotes the beginning velocity, and a = v final t ; andv final = k−k 0 t . By considering v 0 = 0 and substituting these values in Eq. 18, the updating followers slaps can be presented as: (1 + round(rand(0, 1)))] Slime mould algorithm (SMA). The SMA is one of the new swarm meta-heuristic algorithms; it mathematically mimics the propagation wave of slime mould when simulating the best path for connecting foods 39,49 . The details of this method are in 30, and it is proposed in two stages: approaching and warp foods. The details of both stages can be summarized as follows 49 : (a) Approaching food stage: the slime approaching food in this phase based on its odour in the air, and this can be mathematically expressed as follows: where, X represents the position of the slim mould, X A andX B are randomly selected from the mould, X b is the current position related to high odour concentration. t is the current iteration, vc is a parameter gradually decreasing from 1 to zero in a linear form, vc is a parameter ranges from a to − a , where a = arctanh(− t max(t) + 1) . W is the weight of the slim mould and it can be generated and updated based on the fitting accuracy, see 39 . r ∈ (0, 1) , and P = tanh[S(i) − DF], i = 1, 2, 3, . . . , n , here, S(i) is the fitness values of X, and DF is the best fitness over the whole iterations.
(b) Warp food stage: in this stage, the slime behaviour of venous structures can be expressed as: where, LB and UB represent the lower and upper bound of the search range. rand is the random value in between 0 and 1.

Harris hawks optimization (HHO). HHO has been developed by
Heidari et al. 50 as a new optimization technique. It applies a resemblance of Harris hawks cooperative behaviour in optimization solutions 38,50 .
The details of this technique can be found in 50 . In general, HHO depends on three phases that are exploration, transferring, and exploiting 38 . In exploration stage, the position of hawks is determined using the following equation 38 : where, ϒ rand and ϒ prey represent randomly selection hawks and prey's position, respectively. r i denotes a random value in between 0 to 1. ϒ m is the average position.
In transferring phase, the prey energy can be modelled as = 2E 0 (1 − (iter/T)) , where E 0 and Tǫ(−1, 1) , and by determining E, the hawk decides whether to search for exploit the neighbourhood of the solutions. In short, beginning the exploration phase when |E| ≥ 1 , and exploiting the neighbourhood when |E| < 1 . Based on the exploration phase, hawks select soft or hard besiege in applying.

Design algorithms and evaluation
Integration models designs. In the present study, six meta-heuristic OAs, namely PSO, GWO, SOS, SSA, HHO, and SMA are used to optimize the hyper-parameter of LSSVM. These hyper-parameters are γ and σ . Note that the model establishment of LSSVM requires an appropriate setting of its hyper-parameters, including their regularization for constructing an optimum model. The hyper-parameter γ (regularization parameter) and σ (kernel parameter) strongly affect the performance of the LSSVM model, and hence they should be tuned properly for constructing the optimum model in predicting the desired output. In addition to the hyper-parameter, the selection of kernel function also plays an important role. Therefore, proper selection of hyper-parameters of LSSVM at one go is not a trivial task because they should be searched in continuous domains, and hence, there will be an infinite number of parameters sets. Thus, the problem of parameter tuning of LSSVM can be formulated as an optimization problem, and a meta-heuristic OA can solve this.
Considering the above points as a reference, PSO, GWO, SOS, SSA, HHO, and SMA were used to optimize γ and σ of LSSVM and six hybrid LSSVM models LSSVM-PSO, LSSVM-GWO, LSSVM-SOS, LSSVM-SSA, LSSVM-HHO, and LSSVM-SMA, were constructed. The steps of optimizing LSSVM parameters using OAs can be described as follows: (a) initialize LSSVM, (b) set upper and lower bounds of γ and σ , (c) set kernel function, (d) data partitioning, (e) selection of training dataset, (f) initialize OAs, (g) set deterministic parameters of OAs such as, swarm size (n s ), number of iterations (i max ), upper and lower bounds (UB and LB), etc., (h) training of LSSVM algorithm, (i) calculate the fitness function, (j) check and evaluate fitness, (k) obtained optimized values of γ and σ , and (l) testing of hybrid LSSVMs based on obtained values of γ and σ . Figure 1 presents the steps of developing hybrid LSSVM models in the form of a flow chart. Note that, apart from hyper-parameters of LSSVM, the deterministic parameters of OAs also play an important role in hybrid modeling; therefore, they should be tuned appropriately during the course of the optimization process.

Statistical and uncertainty models evaluations.
To evaluate the accuracy of the proposed models, different statistical parameters were applied and evaluated. Correlation statistical parameters, viz., determination coefficient (R 2 ), Nash-Sutcliffe efficiency (NS), and variance account factor (VAF) were calculated to assess the linearity correlation between observed and predicted Mr. Here, when three parameters' values are equal, it means the accuracy of the model is high 51 . Error statistical parameters, viz., root mean square error (RMSE), mean absolute error (MAE), and root mean square error to observation's standard deviation Ratio (RSR), were determined to assess the models' errors. The percentage of error (PE) was calculated to evaluate the model accuracy in predicting Mr. These are the widely used performance indices 10,51-55 , the mathematical expressions of which can be given by:  In addition, Visual evaluation was presented and discussed. In the current study, three visual statistical methods were applied: regression error characteristics (REC) curve, rank analysis, and violin plot. REC curve measures the model accuracy based on the amount of error in the form of squared residuals. The cumulative distribution function of the error between the actual and predicted values is used to determine the model accuracy. The area over the curve (AOC) represents the performance of model accuracy, the model that has a small AOC value is the best model. The rank analysis is another simple visual evaluation technique, it depends on the statistical indices, Eqs. 23-29. In this analysis, the models score a rank from 1 to 6; rank 1 indicates low performance and six refers to the model with good modeling performance. The total rank for the training and testing phases implies which model is the best for modelling the Mr of subgrade soils. Violin plot is another visual analysis technique, it is similar to the box-plot with showing the distribution of probability density of the measured values. This study presented the violin plot for models' errors to compare the best proposed models in the error evaluation. More details for these methods can be found in 55,56 .
For more investigation, the reliability of each proposed models in measuring Mr was assessed through the uncertainty analysis (UA) index. The UA of models can be utilized to test the proposed models under different experimental conditions. Here, for the model error (E i ), the mean of error (MOE) and standard deviation (SD) can be calculated as follows: Moreover, the standard and margin errors, i.e., SE and ME, respectively, are used to calculate the width of the confidence bounds (WCB) 57 . The WCB indicates the upper bound (UB) and lower bound (LB) uncertainty of the proposed models, and it can be determined as: where t denotes the left-tailed inverse of the error distribution. The 95% confidence interval of prediction error can be determined using the values of WCB and MOE. The UB and LB indicate the error range in which approximately 95% of data are located. The lower UA statistical indices, the greater model certainty, indicating a small error in predicting Mr values.

Sensitivity analysis.
To determine the impact of input variables in modelling Mr, the sensitivity of these inputs was conducted. This can provide a guide for using/un-using these variables in the proposed models based on the significance of each variable on the predicted value of Mr. This may help decrease the complexity of proposed models and decrease the measurement cost in future applications. In this study, the cosine amplitude method (CAM) was implemented to assess the impact of input variables 55 . Based on CAM, the strength between Mr and input variables can be determined as follows: where data pairs, x i and x j , of datasets are constructed to measure the strength of the relation. The closer R ij to 100 means more impact corresponding variable has on the Mr value.

Soil sites and datasets
The used dataset of soils in this study was collected from literature presented in the Ohio Department of Transportation 3 . The data were collected from different road construction sites. A total of 891 datasets composed of three types (418 datasets for A-4, 283 datasets for A-6, and 190 datasets for A-7) of subgrade material in Ohio, which are cohesive soils, were used in this study to predict the Mr using the proposed models. The dry side of the optimum, optimum, wet side of optimum, and saturated moister content of soil water conditions were considered in Mr tests, which were performed according to AASHTO standardization 4 . Hanittinan 3 found nine input variables that affect the modeling of Mr, which are percent of soil particles passing through #200 sieve (P200) (fines content), liquid limit (LL), plasticity index (PI), optimum moisture content (OM), soil moisture  Table 1 presents the statistical analysis (maximum (Mx), minimum (Mn), average (M), and standard deviation (SD)) of these variables for the whole datasets, and Fig. 2 demonstrates the histogram and normal distribution of them. From Table 1, it can be seen that the variation between datasets for one variable is high. In addition, from Table 1 and Fig. 2, it can be shown that the distribution of the variables is almost not normal. It can be noted from Table 1 and Fig. 2 a significant variation in the range of dataset for each variable can be detected. This means that a non-linear relationship between inputs and output variables can be detected, and this gives an advantage for using hybrid soft computing techniques in this study compared to traditional techniques.

Results and discussion
Parametric configuration of the developed hybrid models. As stated earlier, the selection of hyperparameters of LSSVM and deterministic parameters of OAs play an important role for constructing the optimum model, therefore, the values of two hyper-parameters ( γ and σ ) were set within a pre-defined wide range of upper and lower bounds. In this study, these parameters' upper and lower bounds were set to (100 and 0.10) and  For instance, the exploration and exploitation constants of PSO were set to 1 and 2, respectively. The value of parameter z in LSSVM-SMA was set to 0.03. It is worth noting that, prior to constructing the models, the main dataset was partitioned into training and testing subsets; among them, the training subset was used to construct the hybrid models, while the testing subset was used to assess the predictive capability of the constructed LSSVM models. The detailed hybrid models constructed with three sets (Set 1, Set 2, and Set 3) are presented in Table 2 for training and training subsets. Herein, the RMSE values of developed models are given in terms of normalized predicted outputs. Right after the model development, they were assessed based on the performance of the testing dataset. It is pertinent to mention that, a model that attained higher prediction accuracy in the testing phase should be accepted with more conviction. As can be seen, all the developed models attained the most accurate prediction when the value of n s was set to 50. Form the results presented in Table 2, it can also be observed that the developed LSSVM-SOS attained the most accurate prediction in the testing phase in all cases, indicating high generalization ability. However, the details of n s , i max , UB and LB of γ and σ , cost function, and the optimum values of γ and σ are presented in Table 3. Herein, the values of these parameters are presented for Set 2 combination of hybrid LSSVM model construction. In addition, the convergence behaviour of three different combinations is presented in Fig. 3. As can be seen, the developed hybrid models converge in less than 20 iterations indicating lower computation cost in all cases. Note that, all the hybrid models were constructed in MATLAB 2015a environment. In the following sub-section, the outcomes of the developed models in predicting Mr of subgrade soils are presented, analyses, and compared. Table 4 presents the statistical evaluation of the proposed models. The correlation between the actual and predicted Mr values is high for all models except LSSVM-SMA model in the training and testing stages. The three statistical correlation parameters (R 2 , NS and VAF) are seen equal for models LSSVM-GWO, LSSVM-SOS, and LSSVM-SSA in the training stage. Furthermore, the model error parameters (RMSE, MAE and RSR) are shown small for the same models in the training (tr) and testing (ts) stages. This means that the accuracy of LSSVM-GWO, LSSVM-SOS, and LSSVM-SSA in modelling Mr is high and can be used as the best models. The percentage error of LSSVM-SOS model is shown low in the training (5.613%) and testing (33) γ and σ = rand × (UB − LB) + LB  Visual and uncertainty models evaluations. Figure 5 illustrates the REC curves for the training and testing stages, and Table 5 presents the AOC of the proposed models. From Fig. 5 and Table 5, it is obviously shown that the accuracy of LSSVM-GWO and LSSVM-SSA is high in the training stage, while the accuracy of LSSVM-SOS in the testing stage is high to predict the Mr value. The worst model in the training and testing stages is obviously the LSSVM-SMA model. The AOC values of LSSV-SMA in the training (0.0113) and testing (0.014) are shown high. The AOC value for the LSSVM-SOS, LSSVM-GWO, and LSSVM-SSA in the testing stage is 0.0016, 0.0017, and 0.0017, respectively. This means that these models' performances are high compared to other models, and the performance of LSSVM-SOS is the best in this study.

Models' performances.
The rank analysis is proposed in Table 6. It can be seen from the table that the rank analyses of LSSVM-PSO From REC and rank evaluations, it can be concluded that the LSSVM-GWO, LSSVM-SOS, and LSSVM-SSA are the best models and can be used in predicting the Mr of subgrade soils. The UA of the developed hybrid LSSVMs was performed and the results are presented in Table 7. To ensure the robustness of the proposed hybrid models, the UA was performed for the testing dataset only. From the information presented in Table 7, it can be seen that the LSSVM-SOS has lower MOE (0.0254), LB (0.215), and UB (0.0294) compared to the proposed models. This means that the accuracy of LSSVM-SOS model in predicting Mr is high at a confidence level of 95%. The SE, ME, and WCB of LSSVM-SOS, LSSVM-GWO and LSSVM-SSA models are the same, this indicates that the three models can be used to estimate the Mr with low uncertainty and high confidence level. A comparison of the whole UA is presented in the rank index; as presented in the table, it can be seen that the LSSVM-SOS has the lowers model errors, with rank 1, among the three models, followed by the proposed models. Therefore, LSSVM-SOS can be used to estimate accurate Mr values of subgrade soils.
Finally, violin plot is presented in Fig. 6 to present the model errors of LSSVM-GWO, LSSVM-SOS, and LSSVM-SSA in the training and testing stages. From Fig. 6, it can be seen that even the LSSVM-SOS model has a maximum error in the training stage, the model error distribution is shown normally. The mean and median of model's errors are the same for the three models. In the testing stage, the maximum errors of models are observed in LSSVM-GWO and LSSVM-SSA. The shape of the violin plot is the same for both models. This means that the performance of both models is the same in predicting Mr. The maximum error of LSSVM-SOS is smaller than www.nature.com/scientificreports/ LSSVM-GWO and LSSVM-SSA, and the variation between the median and mean of the model error is shown smaller than other models. The model's errors distributions are approximately the same for the three models.
These results indicate that LSSVM-SOS is the best to use in predicting Mr of subgrade soils. The developed LSSVM-SOS was compared to some of the previously proposed models designed to detect the Mr of subgrade soils and use the same input variables. Table 8 presents the overall performances of these models in terms R 2 . From Table 8, it can be seen that the new LSSVM-SOS outperforms other historical designed models, and it is a superior alternative to the traditional models in predicting the Mr of subgrade soils. The previous (ANN-GA) and current hybrid models have the same performance, in terms R 2 , accuracy in modelling the Mr with low complexity in the modelling calculations. However, the LSSVM-SOS model outperforms ANN-GA in        Table 7. Results of UA.   Variables impacts on Mr modelling. For better assessment the performance of the developed hybrid models, sensitivity analysis was performed. As stated above, the cosine amplitude method (CAM) 62 was used to perform the sensitivity analysis. Table 9 presents the outcomes of the sensitivity analysis for the different proposed models. In addition, the relative impact of input variables on Mr is presented in Fig. 7 for LSSVM-SOS, -GWO, and -SSA models. It is obviously observed that the impact of all variables in modelling and determining Mr is above 70% (refer to Fig. 7) which means that all variables possess a high impact in modelling Mr values. All variables have approximately the same the impact on the proposed models. Thus, the contribution of all variables cannot be neglected in modelling Mr of subgrade soils. However, the variables fine content (P#200), optimum moisture content, and unconfined compressive strength are shown more significant impact on the Mr, with an impact greater than 80%. These results are in agreement with state-of-the-art studies 58-61 .

Analysis of robustness of LSSVM-SOS.
It is important to note that overfitting is a prevalent problem in data-driven modelling. It means that a data-driven or machine learning model can successfully estimate the desired output during both the training and testing phases, but it can also predict exceedingly odd results for datasets obtained using a completely different design setup. Thus, comparing the overall behaviour of a predictive model to the expected behaviour for a completely different dataset is worthwhile.
In this work, a simulated dataset was constructed to evaluate the robustness, overall behaviour, and expected trend of different input parameters in predicting the Mr of subgrade soils. To generate the simulated datasets, one input parameter was changed while the remaining input parameters remained constant. The details of the simulated datasets are presented in Table 10. Figure 8 shows all of the trends as smooth curves, revealing that when the values of P#200, PI, DS, and σ d increased 63,64 , the Mr of subgrade soils decreased (see Fig. 8a,c,f,i). On the contrary, as the values of LL, OM, SM, q u , and σ 3 increase, the soil Mr increases 63,64 (see Fig. 8b,d,e,g,h). It is worth noting that the LSSVM-SOS model was used to ensure expected trends of the input parameters using a simulated dataset, but real-life analysis may provide different results. Based on the results of the parametric study, the robustness of the proposed LSSVM-SOS model can be established.

Conclusions
In the present study, six swarm intelligence meta-heuristic algorithms (PSO, GWO, SOS, SSA, SMA, and HHO) are applied to optimize the LSSVM parameters for developing a new hybrid technique that can be used for modelling Mr of subgrade soils. 891 datasets of different sites included confining stress, deviator stress, unconfined compressive strength, degree of soil saturation, soil moisture content, plasticity index, percentage of soil particles passing through a #200 sieve, liquid limit, optimum moisture content were used to design and test the proposed models; in addition, the significance of these variables was investigated. The comparison of LSVM-SOS with state-of-the-art models shows that it is extremely efficient in detecting an accurate Mr value for the subgrade soils. However, the future direction of this study may include (a) a comprehensive assessment of the accuracy of the proposed LSSVM-SOS and other hybrid LSSVMs using other datasets of Mr from other regions; (b) an assessment of results of other hybrid models (such as ANNs, ELMs, ANFIS, etc.) constructed with swarm intelligence algorithms 13,65-67 , physics-based algorithms, evolutionary algorithms and human-based OAs; and (c) implementation of different mechanisms such as PSO-based mutation mechanism 67 , adaptive and time-varying acceleration coefficients [68][69][70][71] , Gaussian-based mutation with an exploratory search mechanism 72 , etc. to improve the performance of hybrid models constructed with standard version of OAs. Nonetheless, to the authors' knowledge, this is the first study to apply hybrid LSSVM models created with a specific set of OAs (i.e., swarm intelligence algorithms) to estimate the Mr of subgrade soils.

Data availability
The data used in this study is available in reference 3 and online at https:// etd. ohiol ink. edu/ apexp rod/ rws_ olink/r/ 1501/ 10? clear= 10& p10_ acces sion_ num= osu11 90140 082.  License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.