Antimicrobial data for both B. subtilis and E. coli, used as representative Gram-positive and -negative bacteria, were obtained from the ChEMBL database (version 27) (Gaulton et al., 2012). Using the ChEMBL web server, a dataset was created for each bacteria type by selecting minimum inhibitory concentration (MIC) as the biological/antimicrobial activity measurement. The datasets were then downloaded, and further filtering was performed in the Python environment.

Firstly, compounds with ‘non standard unit for type’ or ‘outside typical range’ in the data validity comments were removed. Then compounds with standard relation values of ‘<’ or ‘>’ were also removed, and duplicates based on compound ‘Molecule ChEMBL ID’ were averaged into one value. At this point the B. subtilis dataset had 9,828 compounds and E. coli had 21,657 compounds, which are hereafter referred to as ‘B-All’ and ‘E-All’, respectively. Since lignin has a chemical composition that only contains carbon (C), hydrogen (H), and oxygen (O) atoms, the datasets were further filtered by keeping compounds with only those atoms. This was performed by searching for compounds with a canonical simplified molecular-input line-entry system (SMILES) string with only C, H, and O atoms (Weininger et al., 1989). The resulting filtering produced a B. subtilis dataset with 768 compounds and an E. coli dataset with 703 compounds, which are hereafter referred to as ‘B-Sorted’ and ‘E-Sorted’, respectively. Finally, to increase the datasets specificity for predicting lignin phenolics, the previously SMILE sorted dataset was filtered for compounds with at least one phenolic ring. This resulted in a B. subtilis dataset with 309 compounds and an E. coli dataset with 278 compounds, which are hereafter referred to as ‘B-Phenolic’ and ‘E-Phenolic’, respectively. Therefore, three datasets for both B. subtilis and E. coli were created with MIC data. Furthermore, MIC values originally determined in µg/ml were converted to micromolar values (µM/ml) and then converted to pMIC (i.e. -logMIC, in molar) for all datasets (Araya-Cloutier et al., 2018). One could consider preprocessing the dataset differently for different ML models, especially for ANN as ANN models can bear more noise than some others. However, the purpose of this study is to compare performance across different ML models using the same dataset. We decided to use the same preprocessing methods across all ML models.

A new dataset of MIC biological activity measurements for lignin monomers against both B. subtilis and E. coli were compiled from published sources. Multidisciplinary databases such as Academic OneFile, Academic Search Complete, EBSCO, and Google Scholar for terms including combinations such as ‘lignin,’ ‘antimicrobial,’ ‘phenolic,’ ‘MIC,’ ‘monomer,’ ‘antibacterial,’ as well as authors with previous work containing appropriate data, were used to find journal articles that contained MIC antimicrobial data for phenolics that can be derived from lignin. In total, 16 compounds were found with MIC data for B. subtilis and 27 compounds for E. coli (listed in Section 3.2). MIC values originally determined in µg/ml were converted to micromolar values (µM/ml) and then converted to pMIC (i.e. -logMIC, in molar) prior to modeling (Araya-Cloutier et al., 2018). The resulting datasets for B. subtilis and E. coli are hereafter referred to as ‘B-Meta’ and ‘E-Meta’, respectively.

The antibacterial activity of 25 lignin derived monomers and three dimers were assessed by monitoring the cell growth (as represented by the optical density at 600 nm, OD₆₀₀) of B. subtilis (NRRL B-354) and E. coli using a spectrophotometry. The full list of compounds and subsequent antimicrobial activity measurements are listed in Section 3.3. The monomers were of analytical quality and purchased from either Sigma Aldrich (St. Louis, MO, USA) or TCI America. The guaiacylglycerol-beta-guaiacyl ether dimer was purchased from TCI America, while 2-(2-methoxyphenoxy)-1-(4-methoxyphenyl)ethanol and 3-hydroxy-2-(2-methoxyphenoxy)-1-(4-methoxyphenyl)-1-propanone dimers were kindly provided by Dr. Mark Crocker at the Center for Applied Energy, University of Kentucky (Song et al., 2018).

Briefly, frozen cultures were first revived in liquid growth media (LB broth, Fisher BioReagents™, BP9723) and allowed to grow at 180 rpm shaking speed for 12 h at 37°C. Afterwards the cells were pelletized, washed, and resuspended in fresh liquid media. To test for the antimicrobial properties, each microbe was cultivated in 96-well plates (run in triplicate) and the OD₆₀₀ was monitored for 24 h with time points taken every 10 min. All wells were brought to an OD₆₀₀ of 0.2 prior to growth, and the phenolics were added to treatment wells to create a final concentration of 1 g/L. To facilitate the solubility of the phenolics in media, all cultures had a final ethanol concentration of 5% (v/v). Two controls were used, one having the 5% ethanol concentration, and one having just microbes and the media. To determine how the phenolics affected microbial growth, the percent change in OD₆₀₀ of the ethanol control during the exponential phase of growth was compared to the growth of the phenolic treatments. This resulted in the percent decrease in growth or Bacterial Load Difference (BLD) for each phenolic treatment (Bouarab-Chibane et al., 2019b), with the formula described in Equation 1:

After obtaining the BLD values for each phenolic, the structures of each compound were converted to canonical SMILES strings using PubChem for use in descriptor calculations. The final experimental datasets for B. subtilis and E. coli are here after referred to as ‘B-Experimental’ and ‘E-Experimental’, respectively.

To calculate the various molecular descriptors, all the compound’s structures in each dataset were converted into canonical SMILES strings (Weininger et al., 1989), if not already provided. These SMILES were then entered into an open-access molecular descriptor calculator software package for Python, RDKit (http://www.rdkit.org). RDKit has a variety of calculatable descriptors that describe a molecule’s lipophilicity (i.e., LogP, LogD), topological indices (i.e., fragment complexity, size, polarity), connectivity indices and different molecular fingerprints (i.e., number of hydroxyl groups, phenolic rings, carboxylic acids etc.). A full list of descriptors and their description is provided in Supplementary Table S1. While it is possible to create models with all the calculatable descriptors, a variety of descriptor selection methods were utilized to improve model accuracy by reducing dimensionality of input space without losing important information.

For the B-All, E-ALL, B-Sorted, E-Sorted, B-Phenolic, and E-Phenolic datasets 200 of RDKit’s available descriptors were calculated. Highly correlated (|r| ≥ 0.8) and constant descriptors were eliminated from the list for each individual dataset. To further reduce the dimensionality of the predictors (descriptors) a principal component analysis (PCA) was performed using scikit-learn (Pedregosa et al., 2011). The number of new principal components to be used was assessed by plotting the number of components vs the percent explained variance, and the number of components that explained 99% of the variance were chosen for each dataset. After the optimal number of principal components were chosen and calculated these values were used as the independent variables for predicting the pMIC values in the subsequent QSAR models. Before modeling, each of the above dataset’s with pMIC and PCA data were randomly split into training (80%) and test (20%) sets three times for cross-validation. We compared and utilized five machine learning algorithms to build the QSAR models for the B-All, E-ALL, B-Sorted, E-Sorted, B-Phenolic, and E-Phenolic datasets. They were the support vector machine (Epsilon-Support Vector Regression), random forest regressor, k-nearest neighbors regressor, decision tree regressor, and neural network regressor (Multi-layer Perceptron regressor) algorithms provided by scikit-learn. The specific settings and parameters used to build each machine learning algorithm are provided below. QSAR models were assessed based on their average coefficient of determination (R²) and root mean squared error (RMSE) based on the predictions made for the three training and test sets.

The best QSAR models constructed from the ChEMBL datasets were further tested for prediction accuracy, by using the meta-analysis datasets as a test set for predicting lignin-specific compounds. Kernel density estimate (KDE) plots using the Seaborn plugin for python were constructed to determine the distribution of each dataset’s pMIC values. Furthermore, the applicability domain (AD) for estimating the reliability in the prediction of new compounds from the ChEMBL datasets were evaluated against the meta-analysis datasets, according to previous work (Golbraikh et al., 2003).

For the B-Meta, E-Meta, B-Experimental, and E-Experimental datasets all 200 of RDKit’s available descriptors were calculated. Pearson’s correlation coefficient (|r| ≥ 0.5) was used to select a fixed subset of predictors (descriptors) best able to predict the antimicrobial activities (either pMIC or BLD) using the ordinary least squares (OLS) regression analysis (Hira and Gillies, 2015). The OLS regressions were performed using Statsmodels (Seabold and Perktold, 2010). As the number of compounds for each of these datasets were very low (less than 30 compounds), the datasets were not separated into training and test sets due to higher risks of chance correlation and overfitting (Araya-Cloutier et al., 2018). For each dataset, the selected descriptors were fed into an OLS regression and backwards elimination was used until the significance of each descriptor coefficient in the model (p-value) was less than 0.05, which identified the best fitting model.

All machine learning models were created using scikit-learn and either the default hyper parameters were used or a number of different parameters through a grid search based exploration of model parameter space was utilized (Lenselink et al., 2017). The final parameters used for the machine learning algorithms that used grid search for QSAR model development are reported in Table 1 .

The support vector machine (SVM) or Epsilon-Support Vector regression is a non-linear regression that calculates an optimal hyper-plane where the distance and error between each data points is minimized (Mei et al., 2005). The SVM performed here used the default parameters provided by scikit-learn. These included a radial basis function kernel, gamma of 1/number of descriptors, parameter cost of 1, and epsilon of 0.1.

Decision tree regressors (DT) are a non-parametric learning method that works by creating a set of binary rules to calculate the target value by dividing the data into subsets that contain data with similar values (Basant et al., 2016). The DT used a grid search to select the optimal maximum depth from 1 to 21 and minimum number of sample leaves from 1 to 100 for each dataset, by fitting the training set and using five cross-fold validations and RMSE to choose the best values. All other parameters utilized the scikit-learn default settings.

The random forest regressor (RF) is an ensemble learning method for non-linear regression analysis, that operates by constructing a multitude of decision trees and outputting the mean prediction of the individual trees (Svetnik et al., 2003). We used all the default parameters provided by scikit-learn, but the number of estimators was increased from the default 100 to 500.

K-nearest neighbor regressions (KNN) are a non-parametric method that stores all available cases and predicts a continuous target based on the similarity measure (distance function) between different features in the same neighborhood (Zheng and Tropsha, 2000). The KNN used a grid search to select the optimal number of neighbors from 2 to 15 for each dataset, by fitting the training set using five cross-fold validations and RMSE to choose the best number of neighbors. The rest of the parameters including the weight function and leaf size utilized scikit-learns default settings.

Neural networks (NN) are brain-inspired algorithms where input features are fed into an input layer, and after a number of nonlinear transformations are performed in a hidden layer, the predictions are generated in an output layer to produce a regression (Lenselink et al., 2017; Camacho et al., 2018; Vamathevan et al., 2019). The ANN relied on most of the default settings provided by scikit-learn’s MLPRegressor neural network. In order to optimize the model’s hyper-parameters, the GridSearchCV function was utilized, with the best parameters being selected based on the best five cross-fold validations and RMSE score. The hyper-parameters chosen to be optimized were the hidden layer sizes, activation type, and learning rate. Specifically, for the hidden layers and number of neurons in each layer was either three layers with 50–100 neurons [(50,50,50), (50,100,50)], or the default setting of a single layer and 100 neurons (100). The activation functions used for the hidden layers were either the rectified linear unit function ‘relu’ or the hyperbolic tan function ‘tanh’. The learning rate schedule for weight updates was either the constant or adaptive learning rates. For further documentation and explanation of the other default settings used to create this model please refer to the scikit-learn version 0.23.2 documentation.

Python (version 3.7.7) was used with the following libraries: RDKit (version 2020.03.6) for the calculation of fingerprints and descriptors, Scikit-learn (version 0.23.2) for all machine learning algorithms and descriptor selection techniques, seaborn (version 0.11.0) with Matplotlib (version 3.3.2) for all figure visualizations, and Pandas (version 1.1.2) for all dataset analysis and manipulation.

The open access database, ChEMBL, was used to develop datasets of compounds with antimicrobial activity measurements (MIC) against both B. subtilis and E. coli. These datasets were used alongside machine learning algorithms to develop QSAR models with the potential to predict the antimicrobial activity of lignin derived phenolics, from compounds that are not lignin specific.

The initial ChEMBL datasets created for B. subtilis (B-All) and E. coli (E-All) contained 9,628 and 21,657 compounds, respectively. These datasets were filtered into two additional subsets, that contained compounds having more similar structures to that of lignin derivatives. The first subset was created by selecting compounds with only C, H, and O atoms, resulting in a B. subtilis dataset with 768 compounds (B-Sort) and an E. coli dataset with 703 compounds (E-Sort). By removing compounds with nitrogenous, chlorine, or fluorine based functional groups, the remaining compounds could have more similar chemical characteristics to that of lignin derivatives. Then those subsets were further filtered by selecting compounds with at least one phenolic ring, resulting in a B. subtilis dataset with 309 compounds (B-Phenolic) and an E. coli dataset with 278 compounds (E-Phenolic). Lignin’s antimicrobial properties have been reported to attribute to its phenolic structures, so it was important to include a subset of compounds that contained only phenolic-based structures. These were the final six datasets used for QSAR model development from the ChEMBL database.

The QSAR models used antimicrobial activities measured in pMIC (-log MIC, in µM/mL) values as the target and the molecular descriptors calculated from RDKit as the variables. Supplementary Table S1 lists the descriptors and molecular fingerprints used in this work. These descriptors represent the lipophilicity (i.e., LogP, LogD), topological indices (i.e., fragment complexity, size, polarity), connectivity indices and functional groups. They are selected based on previous work (Speck-Planche et al., 2012; Svensson et al., 2017; Araya-Cloutier et al., 2018; Bouarab-Chibane et al., 2019b). Supplementary Table S2 (Additional File 2) lists the specific descriptors used for each dataset after pre-processing. Table 1 summarizes the number of descriptors.

Principal component analysis (PCA) was used to reduce the number of descriptors and the dimensionality of the feature space. PCA reconstructs features of a dataset into a new set of uncorrelated features called principal components (PCs). The optimal number of new PCs for each dataset was selected by the number of components that explained 99% of the variance in the dependent variable. Figure 1 shows the number of PCs vs the percent explained variance, and Table 1 summarizes the number that explained 99% of the variance. Since this feature extraction technique creates new independent variables that are less interpretable, the ability to examine how each descriptor influences pMIC is no longer easily obtainable. This is actually beneficial when using the ChEMBL datasets, as we are attempting to predict the antimicrobial properties of lignin with non-lignin based compounds from a data driven perspective, and do not need to understand the exact relationship between these compound’s descriptors and pMIC values. Therefore, the QSAR models were developed from the pMIC and PC values from each dataset using five popular regression-based machine learning algorithms: support vector machine (SVM), random forest (RF), k-nearest neighbor (KNN), decision tree (DT), and neural networks (NN).

The performance summary of five machine learning QSAR models for each ChEMBL dataset is provided in Supplementary Tables S3–S5 , respectively. Each dataset was split randomly into three different training (80%) and test (20%) sets, with each of these same sets being applied to the different model types. The training sets were used to build each machine learning model and the test sets were used for model validation. The metrics used for measuring model performance was the average coefficient of determination (R²) and root mean square error (RMSE) for the three training and test sets. The better performing model is identified as having a high R² and low RMSE value for the average test scores and training scores; thus, Table 2 provides a summary of the best fitting models for each ChEMBL dataset. When comparing models, if one model had a higher R² and lower RMSE for the test sets, but not the training sets, the model with better performance for the test set was chosen, as it is ultimately the more important metric (Bengio et al., 2017). For example, B-All’s best performing QSAR model was the KNN algorithm ( Supplementary Table S3 ), as it had the highest R² of 0.69 for the test set, despite a slightly lower R² for the training sets (0.86) compared to the RF algorithm (0.95). Accordingly, the E-All, B-Sort, E-Sort, B-Phenol, and E-phenol datasets had the most robust QSAR models using the RF, NN, KNN, RF, and KNN algorithms, respectively ( Table 2 ). Not surprisingly the table shows that only NN model has sensible generalization capabilities as there is almost no difference between RMSE for training and test set. The only similar model is E-phenol KNN but it has low prediction quality, so this is most probably effect of high error and low prediction capability of the model in general. The stark difference between training and test group suggests that the model is over-trained or has only interpolation capabilities (like B-ALL KNN and E-ALL RF – huge differences in R²). Future work could test an ensemble system with different techniques to get averaged predications; a possible way to improve the model capability.

A common theme with all the models in each dataset, was that the R² for the test set was always lower than the training set. This could be a sign of model overfitting or unrepresentative data between the training and test sets (Bengio et al., 2017). However, all the models had very low SE values when averaging the R² values of the three different test/train splits for cross-validation, which would suggest compounds are not being underrepresented. The number of independent variables (ICs) used for each dataset were also rather large (80 or 40), which could contribute to overfitting, but they explained 99% of the dependent variable’s variation and when smaller numbers of ICs were used the model’s performance drastically decreased (data not shown). Coupled with the fact that most models used a grid search parametrization technique to fine tune the hyperparameters, these discrepancies may just be a function of the data itself and not with how the models were evaluated or fit. Furthermore, the E-Sort, B-Phenol, and E-Phenol datasets did not have any QSAR models with a R² > 0.6, which is usually needed to describe a truly predictive model (Shahlaei, 2013). Yet, since these datasets are not lignin-specific, the true measure of these model’s performance needs to be evaluated with an additional test set of actual lignin derived compounds.

To this end, the available literature was searched for lignin derived monomers that had reported MIC values against B. subtilis and E. coli. The results from this meta-analysis are reported in Supplementary Table S6 , where 16 compounds were found with MIC data for B. subtilis (B-Meta) and 27 compounds for E. coli (E-Meta). These two datasets were then evaluated as an additional test set for each of the best performing QSAR models found for each ChEMBL dataset, described above. The data is summarized in Figure 2 , where the predicted vs actual pMIC values of the lignin monomers are plotted. It can immediately be seen that none of the ChEMBL QSAR models could accurately predict the lignin monomers. All the models predicted the lignin compounds as having pMIC values roughly less than 2.5, when they are reported as actually having pMIC values greater than 2.5. This suggests these models are grossly underpredicting the pMIC values for the lignin compounds, which would correlate to them having a lower MIC and subsequently greater antimicrobial activity. To understand this, a kernel density estimate (KDE) plot for the ChEMBL and meta-analysis datasets were constructed to visualize the distribution of their pMIC values ( Figure 3 ), and their applicability domains were evaluated ( Supplementary Tables S7, S8 , respectively).

The KDE plots show that the meta-analysis datasets for both E. coli and B. subtilis have pMIC distributions centered around 3–4, while the ChEMBL datasets are centered between 0–2.5. Even though the ChEMBL datasets clearly contain compounds with pMIC values within the distribution of the meta-analysis datasets, they did not lie within the applicability domains (AD) of the ChEMBL datasets. The AD is a useful measure for determining the reliability of a model’s prediction for a new set of compounds. Based on the PCA for each ChEMBL dataset, their ADs were calculated based on the Euclidean distances among all their compounds and a final threshold value is determined (Golbraikh et al., 2003). Then, the same measure is calculated for each of the compounds in the meta-analysis dataset to test if they lie within the threshold of the ChEMBL dataset’s AD. We can see in Supplementary Tables S7, S8 that none of the B-Meta or E-Meta compounds fall within the AD of the ChEMBL datasets. The data show that our data-centric approach to creating datasets of traditional pharmacological compounds that are most similar to potential lignin structures still cannot accurately be used to predict the antimicrobial activity of true lignin derivates. Even though these results are not what the authors had hoped, these data create a more concrete conclusion that a comprehensive dataset of lignin derivatives with antimicrobial measurements needs to be developed. Considering this, QSAR models using actual lignin compounds from the meta-analysis datasets, and an experimentally derived dataset were developed and are discussed in the further sections.

The meta-analysis datasets, used for validating the ChEMBL QSAR models, were applied to develop their own QSARs using ordinary least square (OLS) regressions. Instead of using PCA as a feature extraction technique, univariate feature selection relying on Pearson’s correlation coefficient (r) was employed. Since OLS regressions rely on linear relationships, it made more sense to utilize Pearson’s correlation as it measures the strength of the linear correlation between the independent (descriptors) and dependent variables (pMIC). Therefore, the same 200 molecular descriptors were calculated for the B-Meta and E-Meta datasets, and the descriptors with a r>0.5 were selected to develop the OLS regressions. Subsequently, the selected descriptors were fed into an OLS regression and backwards elimination was used until the significance of each descriptor coefficient in the model (p-value) was less than 0.05, which indicated the best fitting model. It should be noted that this approach is biased for linear relationships and it undercasts models that have non-linear capabilities.

No reliable QSARs using OLS was obtained for the E-Meta dataset (results not shown). This result was not surprising considering the pMIC distribution in E-Meta dataset had three different centers, as shown in the KDE plot ( Figure 3B ). Its variable distribution and small sample size could prevent the QSAR model from capturing any relevant relationships in the feature space (Bouarab-Chibane et al., 2019b). Conversely, even though the B-Meta (16 compounds) dataset was smaller than E-Meta (27 compounds), a more successful QSAR model was developed. To improve model capacity, one could use ANN MLP with two regression outputs – one for E-Meta and the other for B-Meta. Provided the intrinsic relations are not entirely different this would increase dataset and allow training more sophisticated models.

The model for best predicting the antibacterial activity (pMIC) of the lignin monomers in the B-Meta dataset is summarized in Table 3 and Figure 4 . As observed, the selected OLS model showed good predictive power with a R² of 0.759. Three descriptors, SLogP_VSA3, SLogP_VSA5, and fr_AL_COO, were used to develop the best fitting OLS regression model. The SLogP_VSA3 and SLogP_VSA5 descriptors are Molecular Operation Environment (MOE)-type descriptors that bin the output from other descriptor types (i.e., SLogP) and calculate the van der Waals (VDWs) surface area (VSA) of atoms contributing to any specified bin of that output. Thus, SLogP_VSA3 and SLogP_VSA5 calculate the sum of VSA contributions to the lipophilicity measurement SLogP (partition coefficient of compound in two immiscible solvent) within −0.2–0 and 0.1–0.15 bin ranges, respectively. While SLogP and VSA are ‘primary’ descriptors that have a more-or-less interpretable contributions to a compounds mechanism of action, the MOE-type descriptors are intended to be used as model predictors and are not as interpretable (Labute, 2000). Therefore, the negative and positive relationships SLogP_VSA3 and SLogP_VSA5 contribute to the OLS regression, can only be used as a data driven identifier for predicting the pMIC values of lignin compounds. On the other hand, the fr_AL_COO descriptor represents the number of aliphatic carboxylic acid groups in each compound and can directly be used to infer the mechanism of action.

Caffeic, ferulic, sinapic, and p-coumaric acid were the only compounds with an aliphatic carboxylic acid group present in this dataset and they had the lowest observed pMIC values (~3.3). They also represent hydroxycinnamic acid derivatives that are known to have increased antimicrobial properties compared to their more polar hydroxybenzoic acid counterparts (Borges et al., 2013). This is confirmed here by the fact that gallic and protocatechuic acids, with aromatic carboxylic acid groups, had higher MIC values that corresponds to lower antimicrobial activity. Previous work has suggested that hydroxycinnamic acid’s propenoic side chain is responsible for its increased antimicrobial properties, as it facilitates the transport of the molecule through the cell membrane of Gram-positive bacteria (Campos et al., 2003; Nohynek et al., 2006; Borges et al., 2013). Therefore, this explains why an increase in aliphatic carboxylic acid groups correlated to an increase in antimicrobial activity (lower pMIC) for this dataset. Nonetheless, the B-Meta dataset only represents a very small number of lignin monomers and more compounds need to be examined to truly understand or predict the properties that influence their antimicrobial activity.

The antibacterial activity of 25 lignin derived monomers and three relevant dimers were assessed by measuring their BLD or percent inhibition against B. subtilis and E. coli at concentrations of 1 g/L. The BLD values are presented in Table 4 , and they ranged from 2% up to 100%, indicating compounds can be completely inhibitory to both organisms. The 3-hydroxy-2-(2-methoxyphenoxy)-1-(4-methoxyphenyl)-1-propanone lignin dimer was the only compound to show complete inhibition against both B. subtilis and E. coli. Interestingly, the 2-(2-methoxyphenoxy)-1-(4-methoxyphenyl)ethanol lignin dimer only had a BLD value of 66% and 55% for both B. subtilis and E. coli, but its chemical structure differs only by an absence of a methoxy group on β-carbon compared to 3-hydroxy-2-(2-methoxyphenoxy)-1-(4-methoxyphenyl)-1-propanone. Therefore, the presence of this one methoxy group seems to increase the molecules BLD by ~34–45%. Moreover, we can also see from Table 4 that that alkyl chains on the phenolic subunit (4-ethylphenol) and the lignin dimers themselves play an important role in these lignin derivatives antimicrobial properties (i.e., higher BLD values). However, the development of a QSAR model for both organisms will provide an actual statistical relationship between these molecules BLD values and descriptors for more predictive purposes.

The same methods used to the develop the QSAR models for the B-Meta dataset were used for both the B-Experimental and E-Experimental datasets. Where RDKit’s calculated descriptors were chosen based on univariate feature selection (r>0.5) and an OLS regression with backwards elimination was performed until all descriptors had a p-value less than 0.05. The best fitting OLS regressions are summarized in Table 5 and the predicted vs actual BLD values are plotted in Figure 5 . Firstly, the B-Experimental dataset’s OLS regression had greater predictive power with a R² of 0.831 than that of the B-Meta dataset. Four descriptors were used to develop the best fitting OLS regression model: MinABSEStateIndex, PEOE_VSA13, VSA_EState8, and fr_Ar_OH. As stated previously, PEOE_VSA13 and VSA_EState8 are MOE-type descriptors that are intended to be used as model predictors and are not as interpretable for describing the compounds mechanism of action (Labute, 2000). The MinABSEStateIndex is the minimum absolute electrotopological state (E-state) of a skeletal atom, formulated as an intrinsic value plus a perturbation term arising from the electronic interaction and modified by the molecular topological environment of each atom in the molecule (Hall et al., 1991). This descriptor, like the MOE-type descriptors, is used as more of a classification tool for identifying similar compounds instead of describing a feature that could relate to the compounds mode of action. Therefore, while the MinABSEStateIndex, PEOE_VSA13 and VSA_EState8 descriptors show a positive relationship to the lignin compound’s BLD value against B. subtilis. Comparatively, fr_Ar_OH represents the number of aromatic hydroxyl groups in each compound and is better at elucidating their potential antibacterial mechanisms.

The number of aromatic hydroxyl groups can be seen to have a negative relationship with BLD for the B. subtilis data ( Table 5 ). Where gallic acid, catechol, and protocatechuic acid had more than one aromatic hydroxyl group and the lowest BLD values compared to all the other compounds. So, with an increase in the number of aromatic hydroxyl groups there will be a decrease in BLD, correlating to a decrease in the compound’s antibacterial properties against B. subtilis. Bouarab-Chibane et al. (2019a) found a negative relationship between the number of hydrogen donors and the BLD of plant-based polyphenols screened against B. subtilis. Since the number of aromatic hydroxyl groups and the number of hydrogen donors have a direct positive relationship (Bouarab-Chibane et al., 2019a), we can see that in general phenolics with higher overall polarity will have a decrease in antimicrobial properties. This is supported by the experimental data seen here, where highly lipophilic compounds like 4-ethylphenol had high BLD values. However, this model does not provide an explanation for the 3-hydroxy-2-(2-methoxyphenoxy)-1-(4-methoxyphenyl)-1-propanone lignin dimers high BLD value compared to the monomers, highlighting the issue QSAR models can have with limited data sizes and breadth of compound variability.

When examining the OLS regression for E. coli, we can see that the EState_VSA6, VSA_EState3, and fr_aryl_methyl descriptors were used to create the best fitting model. Again, EState_VSA6 and VSA_EState3 are MOE type descriptors that are used for classification-based purposes and cannot be used to infer an influence on the molecules antimicrobial activity. Comparatively, the fr_aryl_methyl descriptor represents the number of aryl methyl groups or an alkyl chain (i.e., ethyl or propyl) on the molecule and it shows a positive correlation with BLD. Thus, as the number of methyl or alkyl chain groups on the phenolic ring increase there is an increase in the BLD or antimicrobial activity of the compounds against E. coli ( Table 5 ). The compounds described by the fr_aryl_methyl descriptor such as 4-ethyl phenol, 4-propyl phenol, and p-creosol are seen to have BLD values of 80–100% for E. coli ( Table 5 ). These compounds are also considered to be more lipophilic given their alkyl groups on the phenolic ring, of which, lipophilicity is already considered an important factor in increasing antibacterial activity against E. coli in the literature (Sikkema et al., 1995; Bouarab-Chibane et al., 2019b). Furthermore, since Gram-negative bacteria have a substantially higher lipid content in their cell wall compared to Gram-positive bacteria (Salton, 1953), it makes sense that the alkylated phenolics were seen to have greater BLD values for E. coli compared to B. subtilis. Additionally, like the B. subtilis OLS regression, E. coli’s model also did not provide inferences as to why the lignin dimer 3-hydroxy-2-(2-methoxyphenoxy)-1-(4-methoxyphenyl)-1-propanone had such a high BLD value.

Overall, when comparing the results from the QSAR models for the meta-analysis and experimental datasets, we can see that the presence of certain compounds and how antimicrobial activity was measured will influence which descriptors play the most important role in describing antimicrobial activity. We saw that the hydroxycinnamic derivatives in the B-Meta dataset drove the negative relationship between the number of aliphatic carboxylic acid groups and pMIC. Comparatively, a higher number of aromatic hydroxyl groups were shown to decrease the BLD in the B-Experimental dataset. Additionally, the E-Meta dataset did not provide an accurate OLS regression model, while E-Experimental’s model showed increasing alkyl groups on the phenolic ring increases BLD values against E. coli. This emphasizes the fact that using different measures of antimicrobial properties and different lignin compounds to develop QSARs for predicting the antimicrobial properties of lignin may lead to different conclusions. It is important to understand the origin of the strains and the cultivation parameters. Meanwhile one could introduce reference compounds across datasets to ensure a more reproducible response. While this is intuitive, the data here provide support for the need in developing a comprehensive and cohesive dataset with lignin derivatives and their antimicrobial properties. Without such a dataset, the ability to accurately predict the antimicrobial potential of lignin and the variety of derivatives that are produced from depolymerization schemes for biorefinery waste stream – lignin valorization is limited.