The considered machine learning approaches can be divided into individual classes. An important model class includes the decision tree based models like Decision Trees (DT), Extra Trees (ET), Random Forests (RF), Gradient Boosting (GB), AdaBoost (ADA), Histogram-Based Gradient Boosting (HGB), Bagging (BAG) and Extreme Gradient Boosting (XGB). In general, decision tree-based models can be seen as non-parametric supervised learning methods which are often used for classification and regression. The value of a target variable is approximated by introducing simple decision rules based on arithmetic mean values for regression approaches as inferred from the data that represent the independent variables. The hierarchy of decision criteria forms different branches in terms of a tree-like structure. The various methods differ in their assumption on their underlying models (Wakjira et al., 2022; Feng et al., 2021). In contrast to a single weak learning model like DT, ensemble methods like ET, RF, GB, XGB, ADA, BAG and HGB consider an ensemble of different weak learning models. The main purpose of ensemble models is the combination of multiple decision trees to improve the overall performance. In more detail, ET and RF are both composed of a large number of decision trees, where the final decision is obtained taking into account the prediction of every tree. In contrast to ET, RF uses bootstrap replicas and optimal split points for decision criteria whereas ET consider the whole original data sample and randomly drawn split points. Overfitting is decreased by randomized feature selection for split selection which reduces the correlation between the individual trees in the ensemble. In addition, in other tree-based ensemble methods, a distinction can also be made between boosting and bagging approaches. Boosting approaches such as GB, XGB and ADA generate a weak prediction model at each step which is added sequentially to the full ensemble model. Such a weighted approach reduces variance and bias which improves the model performance. In contrast, bagging methods such as BAG and HGB generate a weak single DT models in parallel. As follows, Bagging, which stands for boostrap aggregating trains multiple weak learners in parallel and independent of each other. The final individual models are added to the ensemble by a deterministic averaging process which depends on the weights of accurate or inaccurate predictions. The individual models also differ in their definition and consideration of loss or objective functions for predictions in the training phase.

Statistical estimates for the predictive accuracy of the models are usually computed by the root-mean-squared error (RMSE) or the normalized root-mean-squared error (nRMSE) of predictions. The corresponding predicted values Y ̂ n are compared with the actual experimental values Y _n where n denotes the running index in terms of the associated RMSE value R M S E ( Y ̂ , Y ) as defined by R M S E Y ̂ , Y = ∑ n = 1 S Y ̂ n − Y n 2 P (3) with the number of samples P = 188 in our data set. For estimating the model accuracy in comparison with the standard deviation of the target values, one can compute the normalized RMSE values n R M S E ( Y ̂ , Y ) in accordance with n R M S E Y ̂ , Y = R M S E Y ̂ , Y σ Y , Y ̄ (4) with the experimental standard deviation σ ( Y , Y ̄ ) = 1 / P ∑ n P ( Y n − Y ̄ ) 2 and the mean experimental target value Y ̄ = 1 / P ∑ n N p Y n .

SHAP value analysis as developed by Lundberg and Lee (Lundberg and Lee, 2017) is closely related to Shapley values as introduced in game theory (Shapley, 1953). In short, Shapley values provide estimates for the distribution of gains or pay-outs equally among the players (Molnar et al., 2020). Thus, SHAP analysis aims to rationalize a prediction of a specific value by consideration of the feature contributions. The individual values of features in the data set can be interpreted as players in a coalition game (Lundberg and Lee, 2017; Štrumbelj and Kononenko, 2014). In more detail, the algorithm works as follows (Lundberg and Lee, 2017). First, a subset S is randomly selected from all features F. The selected model is then trained on all feature subsets of S ⊆ F. To estimate the feature effects, one model f _S∪{i} is trained with and another model f _S is trained without the feature. The predictions of the two models are compared using the current input f _S∪{i}(x _S∪{i}) − f _S (x _S ), where x _S are the values of the input features in S. Since the effect of withholding a feature depends on other features in the model, the above differences are calculated for all possible subsets S ⊆ F {i}. The Shapley values are then calculated and assigned to the individual features. In more detail, this can be considered as a weighted average of all possible differences with reference to ϕ i = ∑ S ⊆ F i | S | ! | F | − | S | − 1 ! | F | ! f S ∪ i x S ∪ i − f S x S (5) which shows that this approach assigns each feature an importance value that represents the impact of including that feature on the model prediction.

In addition, the Gini feature importance analysis focuses on the mean decrease in impurity of decision-tree based models. For this, one calculates the total decrease in node impurity for a decision tree based model weighted by the probability of reaching that node. The probability of reaching a node can be approximated by the proportion of samples averaged over all trees of the ensemble (Breiman et al., 1984). The Gini index then estimates the probability for a random instance being misclassified when chosen randomly. The higher the value of this coefficient, the higher is the confidence that the particular feature splits the data into distinct groups.

A heatmap of all experimental correlations in terms of the Pearson correlation coefficient r between the target value ln y and all feature values for 188 data points is presented in Figure 1. In addition to diagonal elements, notable correlations in terms of |r| > 0.5 can only be identified for the logarithm of the total protein concentration ln c _t . All other values are smaller than |r| < 0.5 in accordance with negligible correlations. Thus, it can be concluded that cross-correlations between the feature values are of minor importance. The corresponding correlation coefficients for the individual feature correlations with ln y are r = −0.69 (ln c _t ), r = 0.02 (c _D ), r = 0.09 (c _G ), r = 0.18 (q(pH − pI)) and r = 0.28 (c _U ). All values are also visualized in the Supplementary Material. In consequence, non-vanishing positive correlation coefficients for ln y can only be identified for the actual urea concentration and the fraction of protonated titrable groups. In contrast, the total protein concentration c _t shows a strong negative correlation and the correlations for guandidinium hydrochloride and DTT are negligible. Thus, it can be concluded that the total protein concentration dominates the final yield values. The corresponding p-values from a Spearman rank-order correlation coefficient analysis are listed in the Supplementary Material. As can be seen, all values regarding the correlation between ln y and all features are less than p < 0.05 for q(pH-pI), c _U and ln c _t . The corresponding values for c _G and c _D are quite high, which can be explained by the low concentrations, so not many settings can be chosen independently. However, one can conclude that c _G and c _D have a minor effect on the values of ln y. In addition, higher p values between the individual feature correlations are noticeable. This can be understood in terms of the preparation of the dataset obtained from a design of experiment study, which focuses solely on the individual effects of the characteristics on the target values.

As a next step, the corresponding supervised non-optimized machine learning and standard regression methods are assessed in order to predict the corresponding ln y values with regard to a k-fold cross-validation scheme including successive permutations of the training set (Gareth et al., 2013; Wong, 2015). As can be seen in the Supplementary Material, the highest predictive accuracy for a qualitative assessment of the non-optimized models is achieved for gradient boosting and decision tree-based methods. In more detail, histogram gradient boosting (HGB), extra trees (ET), gradient boosting (GB) and random forests (RF) show low nRMSE values between 0.19 and 0.21. The corresponding coefficients of determination between predicted and actual values vary between R ² = 0.72–0.75. Noteworthy, ensemble boosting methods are ideally suited for non-linear regression problems like solubilization processes. Potential reasons for the high accuracy of decision tree-based models were recently published (Grinsztajn et al., 2022). In more detail, decision-tree based models do not overly smooth the solution in terms of predicted target values. Moreover, it was shown that uninformative features do not affect the performance metrics of decision-tree based models as much as for other machine learning approaches. In terms of such findings, it becomes clear that decision-tree based models often outperform kernel-based approaches in terms of predictive accuracies. It has to be noted that ensemble-based methods are also often robust against overfitting issues (Dietterich, 2000). Hence, time-consuming hyperoptimization tunings and validation checks are not of utmost importance.

In contrast to the good predictions of boosting models, standard linear regression methods like least-angle regression (LARS) or Lasso regression (LAS) show a rather poor performance (Supplementary Material). Interestingly, also standard artificial neural networks (ANNs) show a low predictive accuracy. It can be argued that further optimization of hyper parameters may improve the results. All other approaches show a reasonable or even good performance with nRMSE values between 0.26–0.22. As already discussed, the best performance can be observed for advanced decision tree-based ensemble models which reveals the underlying influence of non-linear contributions.

The corresponding predicted and experimental values for ln y from the HGB method are presented in Figure 2. We used a k-fold cross validation approach, where the training data consists of each N-1 data samples with one test data point from the total data set including N samples (Gareth et al., 2013; Wong, 2015). In more detail, each model M _j is trained with the feature data including the samples X = [x ₁, x ₂, …, x _j−1, x _j+1, …, x _N ] and the associated target data Y = [y ₁, y ₂, …, y _j−1, y _j+1, …, y _N ] and predicts the corresponding test target sample Y _j from X _j . As can be seen, X _j and Y _j are not part of the corresponding training data for model M _j . This procedure is repeated for all N models M ₁, M ₂, …, M _j and the corresponding predictions. Notably, a good agreement between the predicted and experimental values for all 0 > ln y > − 1.6 becomes obvious. All values are located within the experimental standard deviation σ(ln y) as denoted by the dashed blue lines. The corresponding good agreement can be rationalized by the large amount of training data in this range. Notable deviations in terms of one outlier with nRMSE > 1 can only be observed for ln y ≤ −1.9 due to missing reference training data. Moreover, the good predictive accuracy is also highlighted by the reasonable value R ² = 0.75 for the coefficient of determination. As shown in the supplementary material, similar conclusions can also be drawn for the prediction of a training data set of 38 samples using an 80/20 ratio split between training and test data. The RMSE values are of comparable quality. In consequence, it can be concluded that solubilization mechanisms and final yield values can be predicted with acceptable accuracy using machine learning approaches.

Due to the reasonable predictive accuracy of certain machine learning models, it can be concluded that the evaluation of feature importances might provide some further insights into the underlying mechanisms of solubilization. In terms of such considerations, we evaluated all training data with the HGB and the ET method in accordance with the SHAP values. The results for the HGB model are shown in the Supplementary Material. As expected, the accuracies of the HGB and the ET model for training data with R ² = 0.91 (HGB), R ² = 0.99 (ET), nRMSE = 0.30 (HGB) and nRMSE = 0.10 (ET) are higher when compared to the predictions from the k-fold cross-validation approach which rationalizes the validity of the following feature analysis.

As a first step, the corresponding Gini feature importance values as calculated from the ET method are presented in Figure 3. As can be seen, the results confirm the dominant influence of the total protein concentration c _t in accordance with the correlation coefficients shown in Figure 1, followed by the urea concentration c _U and the fraction of protonated titrable groups q(pH-pI). With regard to rather low values, the DTT and guanidinium hydrochloride concentrations have negligible influences. It should be noted that the protein’s native structure is characterized by only a very small number of disulfide bonds. This property explains the vanishing influence as well as the relatively low concentration of DTT as a reducing agent in this context. Furthermore, it is well-known that guanidinium hydrochloride is a very potent destabilizing agent. However, in combination with urea it shows a rather complex aggregation behavior around the protein (Oprzeska-Zingrebe and Smiatek, 2021, 2022; Miranda-Quintana and Smiatek, 2021). Accordingly, higher concentrations of guanidinium hydrochloride could probably exert a stronger influence in terms of feature importance. However, it should be noted that this effect is represented here by the somewhat milder denaturation conditions in the presence of urea. In consequence, it can be concluded that the low feature importance of the guanidinium hydrochloride concentration can be rationalized by its corresponding low concentration. These results are also confirmed by the values from the ELI5 analysis which are shown in the Supplementary Material.

Similar conclusions can also be drawn in terms of example decision pathways for an extra trees model with a maximum tree depth of 4 as shown in the Supplementary Material. It becomes obvious that the first decision criterion defines a total protein concentration of ln c _t = 1.39. This value separates between high and low yield branches whereas further criteria based on individual values of q(pH-pI) and c _U are of minor importance and thus lead to different subclassifications.

The corresponding results for the SHAP value analysis with regard to the final yield values ln y from a HGB model are shown in Figure 4. The color-coded beeswarm plot with the associated contributions of SHAP values are depicted in the top panel. The beeswarm plot aims to provide an information-dense summary for the most important features in terms of model predictions for each data instance. Each data instance is represented by a single dot on each feature row. The horizontal position of the dot is determined by the SHAP value of the corresponding feature. In addition, the color illustrates the original value of a feature. Thus, the beeswarm plot highlights the importance or contribution of the features for the whole dataset. As already discussed, it becomes obvious that the largest influence is represented by the total protein concentration, followed by the urea concentration and the amount of protonated titrable groups. Moreover, it can be seen that the total protein concentration has a negative influence on the final yield value. Thus, low values of c _t lead to positive SHAP values and vice versa. Such findings are unique for the total protein concentration while all other factors show a positive correlation. In terms of a molecular understanding, it has to be noted that the solubilization process itself is a rather complex process due to contributions from interfacial phenomena, the composition of the solution as well as further intermolecular mechanisms. Moreover, the individual features and their contributions on the model predictions are ordered in terms of their importance from top to bottom. This ordering is calculated from the mean absolute SHAP value for each feature. In general, such an approach is strongly determined by the broad average impact of the feature while rare maximum or minimum values do not contribute significantly. The corresponding mean absolute SHAP values are presented in the bottom of Figure 4 for reasons of consistency. It clearly can be seen that the total protein concentration dominates the feature importance, followed by the urea concentration and q(pH-pI). The contributions from c _G and c _D are of minor importance in agreement with the p values from the previous correlation analysis.

For a more detailed analysis, we present the corresponding SHAP value dependency plots for the total protein concentration in Figure 5. In general, SHAP dependency plots as shown in Figures 5, 6 highlight the effect of a single feature on the model predictions. Each dot corresponds to a single prediction from the dataset and the position on the horizontal axis denotes the corresponding actual value of the feature. In contrast, the vertical axis shows the SHAP value for that feature and its impact on the prediction. In general, a slope along the points in the dependency plots enables to identify positive, negative or no correlations with the corresponding target parameter. It clearly can be seen in Figure 5 that the aforementioned negative correlation (Figure 4) can be interpreted as an inverse relation between the total protein concentration and the final yield value ln y. Hence, for increasing total protein concentrations, one can observe a linear decrease of the SHAP values. The corresponding mechanism can be explained as follows. In general, the yield is defined by y = N f N t (6) which corresponds to the ratio between the number of free (solubilized) protein chains N _f and the total number of chains N _t . With regard to the fact that inclusion bodies show a rather poor solubility, one can assume that individual inclusion bodies aggregate in order to form larger moieties. Hence, the dissolution of free chains mainly occurs at the solvent accessible surface area of these compounds in agreement with previous assumptions (Walther et al., 2013).

The number of free protein chains can be written as N f = c p R I 2 d with the local concentration of proteins in the aggregated inclusion body c _p , the corresponding spherical radius R _I of the aggregate and the penetration or dissolution depth d. Moreover, it is assumed that the inner region of the compound remains unaffected by dissolution such that N t ≈ c p R I 3 with R _I ≫ d. In consequence, we can rewrite Eq. 6 according to y ∼ R I 2 d R I 3 ∼ d R I (7) which reveals that aggregated inclusion bodies with larger radii result in lower yield values. Furthermore, it is assumed that d is constant after a fixed time interval. With regard to the fact that the radius scales as R I = N t 1 / 3 / c p 1 / 3 , one obtains after insertion into Eq. 7 the following relation y ∼ d R _ I ∼ d c p 1 / 3 N t 1 / 3 (8) which clearly shows that the final yield value depends inversely on the total protein concentration N _t ∼ c _t as represented by ln y ∼ − ln(c _t ) in agreement with Figure 5. Here, we assume that d is constant after a fixed time interval and that aggregates show a highest packing fraction leading to fixed values of c _p .

Moreover, it is well known that the presence of urea induces a structural destabilization of proteins. Recent articles rationalized this phenomena with preferential binding and exclusion mechanisms (Smiatek, 2017; Oprzeska-Zingrebe and Smiatek, 2018; Miranda-Quintana and Smiatek, 2021). As can be seen in Figure 6, one observes increasing SHAP values and thus a positive trend of ln y for increasing urea concentrations. Moreover, it can be seen that for urea concentrations larger than 8 mol/L, a saturation behavior becomes evident. As an explanation, we refer to co-solute induced destabilization effects as discussed in Oprzeska-Zingrebe and Smiatek (2018); Smiatek (2017). In more detail, it is assumed that the ratio of destabilized and stable proteins K _cs can be written as K c s = K 0 ⁡ exp a 33 Δ ν 23 (9) with the derivative of the thermodynamic activity a ₃₃, the difference in the preferential binding coefficients Δν ₂₃ and the ratio of destabilized and stable proteins K ₀ in absence of any co-solutes (Smiatek, 2017; Oprzeska-Zingrebe and Smiatek, 2018; Smiatek et al., 2018). For certain proteins, it was discussed that the partial molar volumes and the solvent-accessible surface area upon unfolding do not change significantly according to Δν ₂₃ = c _U ΔG ₂₃ (Smiatek et al., 2018; Krishnamoorthy et al., 2018a), such that the previous relation can be approximated by K c s ≈ K 0 ⁡ exp a 33 c U Δ G 23 (10) with the difference in the Kirkwood-Buff integrals ΔG ₂₃ (Oprzeska-Zingrebe and Smiatek, 2018). In addition to stable and destabilized proteins, one can apply the same relation for the fraction of dissolved and bound proteins (Smiatek et al., 2018). Under the assumption that N _f ≪ N _t , it thus follows that K _cs ≈ y and K ₀ ≈ y ₀. For high co-solute concentrations, it was further discussed that a ₃₃ → 0 due to stability conditions (Krishnamoorthy et al., 2018b). Hence, the exponential factor in Eq. 10 can be linearized according to y ≈ y 0 1 + a 33 c U Δ G 23 (11) which highlights the linear contribution y ∼ c _U between the urea concentration and the final yield value in good agreement with Figure 6.

In our previous discussion, it was also highlighted that the fraction of protonated titrable groups has a positive influence on the final yield values. In accordance with Eq. 2, it becomes clear that q(pH-pI) decreases with increasing pH values. Moreover, it is known that the isoelectric point with pI = 8.4 corresponds to a net-uncharged protein. The clear distinction between pH values below and above the pI value in terms of the SHAP dependency plot can be seen in Figure 7. In more detail, q(pH − pI) < 0.5 as represented by pH < 8.4 results in negative SHAP values and thus lower yields and vice versa. The reason for this observation might be related to electrostatic repulsion between the charged groups of the protein (Smiatek et al., 2020). Thus, we assume that for high or low pH values, depending on the net amount of basic or acidic groups in the protein, the electrostatic repulsion fosters the dissolution of the protein in order to minimize non-favorable interactions. As also illustrated in Figure 7, significant contributions for a net-uncharged protein at pH = 8.4 leading to q(pH − pI) = 0.5 are absent.

The combination of the previous considerations in terms of Eqs 2, 8, 11 results in y ∼ y 0 1 + a 33 c U Δ G 23 q p H − p I c t 1 / 3 (12) which condenses all previous results in one analytic expression. In order to assess the validity of Eq. 12, we plotted all experimental data points onto a master curve as shown in Figure 8. In order to minimize fluctuating electrostatic repulsions, we chose nearly constant and high pH values with pH ≥ 10 as well as high urea concentrations with c _U ≥ 7.55 mol/L. The corresponding 27 data points are then scaled in terms of a normal distribution and compared to Eq. 12. As can be seen in Figure 8, the theoretically predicted values nicely follow the proposed scaling relation. It has to be noted that the corresponding influences for different proteins and modalities upon solubilization may vary, such that the obtained results are not generally transferable without proper assessment. Nevertheless, we have proven strong evidence that explainable machine learning approaches provide deeper insights into the molecular mechanisms and correlations of solubilization processes.