Computational Thermodynamic Analysis of the Interaction between Coagulants and Monosaccharides as a Tool to Quantify the Fouling Potential Reduction in the Bioﬁlm Membrane Bioreactor

: The membrane bioreactor (MBR) and the bioﬁlm membrane bioreactor (BF-MBR) are among key solutions to water scarcity; however, membrane fouling is the major bottleneck for any expansion of these technologies. Prepolymerized aluminum coagulants tend to exhibit the greatest extent of fouling alleviation, with the reduction of soluble microbial products (SMPs) being among the governing mechanisms, which, nevertheless, has been poorly understood. This current study demonstrates that the investigation of the chemical coordination of monosaccharides, which are the major foulants in MBR and BF-MBR, to the main hydrolysis species of the prepolymerized aluminum coagulant, is among the key approaches to the comprehension of the fouling mitigation mechanisms in BF-MBR. Quantum chemical and thermodynamic calculations, together with the multivariate chemometric analysis, allowed the team to determine the principal mechanisms of the SMPs removal, understand the thermodynamic patterns of fouling mitigation, develop the model for the prediction of the fouling mitigation based on the thermodynamic stability of the inorganic-organic complexes, and classify these complexes into thermodynamically stable and less stable species. The results of the study are practically signiﬁcant for the development of plant surveillance and automated process control with regard to MBR and BF-MBR systems.


Introduction
Membrane bioreactor (MBR) and biofilm membrane bioreactor (BF-MBR) are advanced solutions for the problem of water scarcity, which have been recognized as highly competitive technologies when applied in water reuse schemes [1][2][3][4]. However, membrane fouling remains the major barrier for any MBR and BF-MBR expansion [5,6].
With regard to membrane fouling mechanisms, cake and gel layer formation and membrane pore blockage were identified as major contributors to any filtration resistance [7]. The formation of the cake layer is mainly attributed to the deposition of suspended solids, whose size is bigger than the membrane pores, onto the membrane surface or onto the sealed pores, with a subsequent stacking [7,8]. The gel layer is the matrix, which consists of highly concentrated solutes and macromolecular species, deposited at the membrane surface. the thermodynamic properties of the Keggin Al 13 sulfate and selenate molecules in the crystal state, in order to identify the link between their structure and their reactivity. The highly exothermic enthalpies of formation of the Keggin Al 13 clusters from the elements and oxides, and the enthalpies of the relevant 5N solutions, were reported.
The interaction between the aluminum hydrolysis species and glucose was studied by He et al. [44], who focused on the investigation of the coordination of β-D-glucopyranose to [Al(OH)(H 2 O) n ] 2+ and [Al(OH) 2 (H 2 O) n ] + ions through the Density-functional method. The formation of the double-O-ligand coordination complex was found to be thermodynamically favorable if it was formed through two O 4 -O 6 adjacent oxygen atoms in the β-D-glucopyranose moiety.
However, no computational thermodynamic studies can be found on the chemical coordination of the monosaccharides to the polymeric hydrolysis species of the prepolymerized aluminum chloride, especially Al 13 n+ , which constitute the majority of its hydrolysis species [45,46] and have the highest coagulation potential [47], thus being primarily responsible for the flux enhancing ability of the relevant coagulant. The current work aims to fill the existing conspicuous gap in the scientific knowledge by developing a strategy, which combines the quantum chemical model and thermodynamic calculations with the multivariate chemometric approach to identify the thermodynamically favorable pathways of the formation of inorganic-organic complexes, and the thermodynamic patterns of fouling mitigation during the application of the prepolymerized aluminum coagulant of the medium basicity in BF-MBR.

Materials and Methods
The calculations of the thermodynamic parameters of the formation of the individual reactants and the relevant complexes were carried out, applying the semi-empirical computational molecular orbital method-Parameterized Model number 3 (PM3) using HyperChem TM 8.0.6 software (Hypercube Inc., Waterloo, ON, Canada). PM3 is based on the neglect of diatomic differential overlap approximation. This parametrization procedure enables the acquisition of fully optimized molecular geometries and the calculation of the heats of formation, dipole moments and ionization potentials. In comparison with the MNDO (Modified Neglect of Diatomic Overlap) and AM1 model (Austin Model 1), the PM3 method is characterized by lower average absolute errors when calculating standard enthalpies, and provides more precise calculations [48]. The PM3 method provides the enthalpies of formation with the weighted total mean absolute deviation ±10.9 kJ/mol [49].
In order to simulate the interactions between the monosaccharides and Al 13 -complex, the relevant molecules were built and solvated using 216 water molecules and the periodic solvent box of the size 18.7 × 18.7 × 18.7 Å for the monosaccharides; and the periodic box of 31.3 × 31.3 × 31.3 Å using 1010 water molecules for Al 13 and Al 13 -monosaccharide complexes. The minimum distance between the solvent and the solute atoms was kept at 2.3 Å. Full geometry optimization using a Polak-Ribiere algorithm was performed for every analyzed compound, in order to reach the conformation of the lowest potential energy. The restricted Hartree-Fock method was applied to calculate the spin interactions in every compound. The standard temperature was set in all of the simulations (298.15 K). The results of the quantum chemical calculations are reported with the accuracy of the convergence parameter, i.e., the self-consistent field energy (SCF), equal to ±0.0418 kJ/mol.
The calculation of the main thermodynamic parameters of the reactions (i.e., standard enthalpy change (∆H o ), standard entropy change (∆S o ), and standard Gibbs energy change (∆G o )) between the selected monosaccharides and Al 13 complex, resulting in the formation of the relevant Al 13 -monosaccharide complexes, was performed according to Hess's law (Equation (1)), its extension to entropy (Equation (2)) and Gibbs energy (Equation (3)), and the Gibbs-Helmholtz equation (Equation (4)).
where ∆H o f ,i is the standard enthalpy change of formation; S o f ,i is the standard entropy of formation; ∆G o f ,i is the standard Gibbs energy change of formation of the individual reactants and products. The statistical investigation of the calculated thermodynamic parameters of the studied systems was conducted using PLS regression analysis and cluster analysis. The statistical software, The Unscrambler®X10.3 (CAMO Software AS, Oslo, Norway), was used for this purpose.

The Main Identified Foulants
Monosaccharides (C n H 2n O n , n = 3-6) and their derivatives are the main building block molecules of polysaccharides, which determine the characteristics of the latter ones. Therefore, it was decided to investigate the behavior of the monosaccharides in terms of their interaction with the selected coagulant.
Different studies investigate the presence of the monosaccharide species in the mixed liquor and the fouling layers of the membrane bioreactor (MBR) systems, as well as in the biofilm matrixes (Table 1).
Notes: 1 if the OH group located at the bottom-most asymmetric centre (the carbon atom second from the bottom) is on the right side in the Fischer projection, the monosaccharide (or its derivative) belongs to D-sugars, if the OH group at the bottom-most asymmetric centre is on the left side in the Fischer projection, the compound is L-sugar [56]; 2 "+" stands for the presence of the selected monosaccharide in the relevant source.
Water 2019, 11, 1275 6 of 30 According to Table 1, such monosaccharides as glucose, galactose, glucuronic acid, glucosamine, rhamnose, fucose, arabinose, and mannose (marked with the blue rectangle) are the most typical constituents of the polysaccharides in the mixed liquor, gel, and cake layers in the MBR systems. Since the vast majority of them were also identified in the biofilm matrices [52], in the present study, it is assumed, that the highlighted monosaccharides make up the polysaccharides in the BF-MBR system.
Taking into account the above-mentioned references, a couple of nuances can be singled out. First, the relative proportions of the monosaccharide species in the structures of the polysaccharides in the gel and cake layers of the MBR system tend to vary, depending on the applied solids retention times (SRTs), which is demonstrated by Silva et al. [50]. Besides, the authors point out that the variety of monosaccharides increases with the decrease of SRT and the intensification of biofouling (all of the dominant ones at different SRTs monosaccharides lie in the blue rectangle in Table 1). Second, as reported by Miyoshi et al. [55], certain monosaccharides might have a higher affinity to the membrane surface, and can thus cause more severe irreversible fouling than others. The clear indicator was the difference in the monosaccharide content in the fouling layer and in the mixed liquor. According to the results [55], glucose, galactose, rhamnose and mannose were found to be dominant in the matrix of the fouling layer.
Monosaccharides in the cyclic form are characterized by an active center, which is called the glycosidic hydroxyl group. The higher reactivity of the glycosidic hydroxyl group can be explained by the influence of the ether-type oxygen atom (between C 1 and C 5 ) (Figure 1a), which partially shifts the electrons from the contiguous C-O bond (in the C 1 position) to its own orbitals.
Hence, the shift of electron density increases the polarity between the C 1 carbon and the -OH group, making this hydroxyl group more chemically active [68]. The shift of the electrons from the carbon atom in position 1 (atom C 1 ) to the ether-type oxygen atom in the structure of β-D-glucopyranose, is clearly indicated by the areas of the relatively higher electrostatic potential (around 1.25 e/a 0 ) and the increased net positive charge of C 1 (+0.194) in comparison to the other present carbon atoms, which is demonstrated in Figure 1b. Figure 1b demonstrates the computationally generated electrostatic potential map of the β-D-glucopyranose molecule, which visualizes the charge/electron density distribution within the molecule, represented as the isosurface. According to the gradation of color in the electrostatic potential map, the asymmetrical distribution of the electron charge density is clearly indicated in the O-H groups. The hydrogen atoms have a low electron density, and hence a net positive charge and positive electrostatic potential, as shown by the bright green color, while the oxygen atoms have a high electron density, and thus a net negative charge and negative electrostatic potential, indicated by the deep purple color. Meanwhile, the C-O groups are characterized by the accumulation of the net positive charge and positive electrostatic potential at the carbon atoms and the negative charge and negative electrostatic potential at the oxygen atoms. The accumulation of the net negative charge at the oxygen atoms can be explained by the fact that oxygen is the element of the second highest electronegativity (χ) in the periodic table (χ O = 3.5, vs. χ C = 2.5 vs. χ H = 2.1 on the Pauling scale), and therefore its atoms have the highest relative ability (after fluorine) to attract the electrons of other atoms to which it is bonded [69].
The electrostatic potential map indicates a significant part of the intermediary electrostatic potential regions, whose colors are not completely green or completely purple. This effect is explained by the nature of observed bonds in the structure of β-D-glucopyranose, as well as in the structure of the other selected monosaccharides, which is covalent polar. Meanwhile, the polarity of the O-H bonds is higher than the polarity of the C-O bonds-∆χ OH = 1.4 > ∆χ CO = 1.0.
The electrostatic potential surfaces of the other investigated monosaccharides have the character similar to the one of β-D-glucopyranose and are demonstrated in Figure S1. The β-D-glucopyranose molecule was selected for the description since it is the most stable aldohexose, and is among the most abundant monosaccharides in MBR [50,56].
The total dipole moment of the β-D-glucopyranose molecule is 1.639 D, with the highest contribution from the X and Z vector components (Figure 1a). Monosaccharides in the cyclic form are characterized by an active center, which is called the glycosidic hydroxyl group. The higher reactivity of the glycosidic hydroxyl group can be explained by the influence of the ether-type oxygen atom (between C1 and C5) (Figure 1a), which partially shifts the electrons from the contiguous С-O bond (in the C1 position) to its own orbitals.
Hence, the shift of electron density increases the polarity between the C1 carbon and the -OH group, making this hydroxyl group more chemically active [68]. The shift of the electrons from the carbon atom in position 1 (atom C1) to the ether-type oxygen atom in the structure of β-Dglucopyranose, is clearly indicated by the areas of the relatively higher electrostatic potential (around 1.25 e/a0) and the increased net positive charge of C1 (+0.194) in comparison to the other present carbon atoms, which is demonstrated in Figure 1b. The glycosidic hydroxyl group located at the C1 atom (marked with the red oval); (b) electrostatic potential (e/a0), mapped onto an isosurface of the total electron density of 0.13 e/a0 3 (Notes: the numerator stands for the electron charge "e" (1.6022 × 10 −19 coulombs), and the denominator is the first Bohr radius "a0" (5.62918 × 10 −11 m).). Figure 1b demonstrates the computationally generated electrostatic potential map of the β-Dglucopyranose molecule, which visualizes the charge/electron density distribution within the molecule, represented as the isosurface. According to the gradation of color in the electrostatic potential map, the asymmetrical distribution of the electron charge density is clearly indicated in the O-H groups. The hydrogen atoms have a low electron density, and hence a net positive charge and positive electrostatic potential, as shown by the bright green color, while the oxygen atoms have a high electron density, and thus a net negative charge and negative electrostatic potential, indicated by the deep purple color. Meanwhile, the C-O groups are characterized by the accumulation of the net positive charge and positive electrostatic potential at the carbon atoms and the negative charge and negative electrostatic potential at the oxygen atoms. The accumulation of the net negative charge at the oxygen atoms can be explained by the fact that oxygen is the element of the second highest electronegativity (χ) in the periodic table (χ O = 3.5, vs. χ C = 2.5 vs. χ H = 2.1 on the Pauling scale), and therefore its atoms have the highest relative ability (after fluorine) to attract the electrons of other atoms to which it is bonded [69].
The electrostatic potential map indicates a significant part of the intermediary electrostatic potential regions, whose colors are not completely green or completely purple. This effect is explained by the nature of observed bonds in the structure of β-D-glucopyranose, as well as in the structure of the other selected monosaccharides, which is covalent polar. Meanwhile, the polarity of the O-H bonds is higher than the polarity of the C-O bonds-Δχ OH = 1.4 > Δχ CO = 1.0.

Dominant Hydrolysis Species with Regard to the Medium Basicity Prepolymerized Aluminum Chloride
Prepolymerized aluminum coagulants are characterized by the dominance of the species with the charges higher than the charges of the majority of the species of non-prepolymerized inorganic coagulants. This feature enhances the coagulating ability of the prepolymerized aluminum coagulants and simplifies the operation process, thus making them more advantageous than their non-prepolymerized counterparts [35]. One of the main factors which influences the dominance of certain polyaluminum hydrolysis species is the basicity of the prepolymerized aluminum coagulant [70]. According to the previous work by Kulesha et al. [31], prepolymerized aluminum coagulant with the medium basicity (OH/Al 1.3) exhibits the best potential to mitigate fouling in the BF-MBR system at the optimum pH range 5.5-6 (acidic), which was mainly attributed to the highest bearing charge concentration among the studied species.
Al 13 7+ ( Figure S2) is considered the most stable aluminum complex in the prepolymerized aluminum chloride (PACl) solution [71,72], which is mainly attributed to a surface positive charge with the π-electrons delocalized in the six-member (hexagon-like) ring structures [73]. Meanwhile, the surface charge decrease makes the subsequently formed Al 13 7 − n species less stable [74]. The presence of Al 13 7+ complexes in polyaluminum chloride solutions is the prime contributor to their  12 ] 7+ ) in the medium basicity coagulant of OH/Al 1.5 at pH ≈ 3.8 at the concentration 0.1 M Al, was reported by Bottero et al. [77], during the analysis of the nuclear magnetic resonance (NMR) spectra (~71% out of the total aluminum concentration) and computational analyses based on Glueckauf's formula and the Debye-Hückel law (~80-88% out of the total aluminum concentration). The experimental potentiometric titration analysis and the relevant model based on its results indicated the presence of the [Al 13 O 4 (OH) 24 (H 2 O) 12 ] 7+ complex, which consisted of the symmetrical tetrahedrally-coordinated aluminum ion at the center of the structure with the weak/non-existent electric field gradient, surrounded by twelve octahedral aluminum ions with the relatively high electric field gradient and the potentially distorted structure.
Feng et al. [78], who studied the speciation of different prepolymerized aluminum coagulants by applying electrospray ionization mass spectrometry, reported that the Al 13 3+ species was one of the main components of these prepolymerized aluminum coagulants, which was assumed to be directly transformed from Keggin-[Al 13 O 4 (OH) 24  Pernitsky and Edzwald [70], based on the experimental solubility data for the sulfated prepolymerized aluminum chloride coagulant of the medium basicity, hypothesize either the presence of the combination of monomers (Al(OH) 2 + , Al(OH) 2+ ) and Al 13 7+ , or the occurrence of some other aluminum species in the studied system at pH 5-7.
In the study by Rämö et al. [79], who investigated the distribution of the polyaluminum species in the 1 mum Al solutions of prepolymerized aluminum chloride with the medium basicity (OH/Al 1.3), based on the mass spectrometry (MS) analysis, the dominance of Al 13 2+ and Al 13 3+ compounds, which made up 72% out of the total ion count at pH 4.7, is reported. No detection of [Al 13 O 4 (OH) 24 (H 2 O) 12 ] 7+ could be explained by any availability of the counter anions in the system and the drying droplet of the specimen. Some larger formations, for example, Al 26 , were also observed; however, the Al 13 complex is the dominant species in the system. On the contrary, the monomers, dimers, and trimers of aluminum were completely absent [46]. According to mass spectrometric studies, reported by Sarpola [45], the PACl with the medium basicity in the pH range 4.65-6.46 is characterized by the following protonated open form of Al 13 -[Al 13 O 4 (OH) 29 ] 2+ . Due to the high stability of the [Al 13 O 4 (OH) 24 (H 2 O) 12 ] 7+ complex, it is suggested that the complex in its usual form is less likely to participate in any reactions. Meanwhile, it is hypothesized [46] that after the formation of [Al 13 O 4 (OH) 24 (H 2 O) 12 ] 7+ in the aqueous solutions, the four-coordinate oxygen atoms in each ring, which are shared with the central aluminate ( Figure S2), get "over-bonded", and thus "loose" the bond with one of the other twelve aluminum atoms that they are connected to, getting protonated and exposing the open chains or ring substructures to water. Hence, four rings with the asymmetrical charge division are produced, giving rise to the active centers, which can potentially attract the negatively charged systems and react with them [46].
Seichter et al. [80], who investigated the species of prepolymerized aluminum chloride, formed by hydrolysis and condensation, assigns the detected complex cations to the octahedral structure [Al 13 (OH) 24 24 ] 15+ , which is identified as the other important tridecameric cation, in addition to the Keggin (tetrahedral) type. This complex, as the complex [Al 13 O 4 (OH) 29 ] 2+ introduced by Sarpola [45], has the planar core, and thus is potentially much more likely to participate in reactions with the organic matter, as well as undergo further polymerization. However, the formation of this polycation structure was expected to occur autonomously from the Keggin-type cation.
Based on the discoveries by Sarpola [45] and all above-mentioned findings, the present study assumes that during coagulation, applying the prepolymerized aluminum chloride with medium The planar open structure of the Al 13 3+ complex was introduced by Sarpola [45], and taken as the basis for building Al 13 2+ . Meanwhile, the geometry optimization analysis conducted in the present work indicates that the minimum of potential energy can only be reached if the molecule has a non-planar conformation, as represented in Figure 2.
The pairs of aluminum atoms in the Al13 2+ complex, which can potentially participate in the interaction with the listed active centers of the monosaccharides, are specified in Section 3.4.
It was assumed that the relevant active centers of the monosaccharides first got deprotonated, i.e., ionized, and then they participated in the chemical coordination with the Al13 complex. As a result of these interactions, two water molecules are dehydrated, which was taken into consideration during the relevant calculations. The enthalpy and Gibbs energy of formation of the liquid water molecule specified by Dean [81] were used in this study. The entropy change for water molecules was calculated from the reference data [81] according to Equation (2).
The assessment of the spontaneity of the potential interactions between the selected monosaccharides and Al13 2+ complex was performed through two steps: First, a series of quantum chemical simulations were conducted, focused on the determination of the thermodynamic parameters of the formation of the individual reactants and the relevant complexes-standard According to the charge/electron density distribution within the molecule, represented as the isosurface of the Al 13 -complex (Figure 2a), aluminum atoms have a low electron density, a net positive charge and a positive electrostatic potential, as shown by the bright green color; meanwhile, the oxygen atoms have a high electron density, and thus a net negative charge and a negative electrostatic potential, indicated by the bright purple color. According to the represented electrostatic potential map, all aluminum atoms, except for the central tetrameric aluminum, are the potential active centers during the interaction with the foulants. Meanwhile, based on the charge balance calculations, it was identified that Al no. 64 ( Figure 2b) has a slight net negative charge (−0.023), and that Al no. 67 has almost no charge (0.004); thus, both are less likely to participate in these reactions. Consequently, the aluminum atoms, highlighted with the light green color in Figure 2b, are the potential active centers of the Al 13 -complex. The total dipole moment of the Al 13 complex is 7.25 D. The analysis of the components of the total dipole moment demonstrates the maximum contribution from the Z vector component (Figure 2b).
The nature of the observed bonds in the structure of the Al 13 -complex is ionic, since the polarity of the O-Al bonds (∆χ AlO ) is equal to 2 (χ O = 3.5, χ Al = 1.5 on the Pauling scale). The great difference in electronegativities is also indicated by the significant differences in the electron density at the aluminum and oxygen atoms, which is demonstrated by the entirely green and entirely purple regions, respectively, with no regions of intermediary electrostatic potential.
The pairs of aluminum atoms in the Al 13 2+ complex, which can potentially participate in the interaction with the listed active centers of the monosaccharides, are specified in Section 3.4. It was assumed that the relevant active centers of the monosaccharides first got deprotonated, i.e., ionized, and then they participated in the chemical coordination with the Al 13 complex. As a result of these interactions, two water molecules are dehydrated, which was taken into consideration during the relevant calculations. The enthalpy and Gibbs energy of formation of the liquid water molecule specified by Dean [81] were used in this study. The entropy change for water molecules was calculated from the reference data [81] according to Equation (2).
The assessment of the spontaneity of the potential interactions between the selected monosaccharides and The standard Gibbs energy change of the process of formation of Al 13 -monosaccharide complex (interaction between the selected monosaccharides and Al 13 2+ complex), which is the main indicator of the spontaneity in chemical reactions, is represented in Figure 3. The analysis of the acquired charts ( Figure 3) demonstrates that all of the processes which result in the formation of the suggested Al 13 -monosaccharide complexes have a negative standard Gibbs energy, which indicates that the processes of their formation are spontaneous, and hence thermodynamically favorable. Meanwhile, the thermodynamic stability of the formed Al 13 -monosaccharide complex highly depends upon the active centers of the monosaccharides and the Al 13 2+ complex, which participate in the chemical coordination process, and the nature of the monosaccharide. The following complexes are identified as the most thermodynamically stable, since they exhibit the highest negative values The standard Gibbs energy change of the process of formation of Al13-monosaccharide complex (interaction between the selected monosaccharides and Al13 2+ complex), which is the main indicator of the spontaneity in chemical reactions, is represented in Figure 3.   The analysis of the acquired charts ( Figure 3) demonstrates that all of the processes which result in the formation of the suggested Al13-monosaccharide complexes have a negative standard Gibbs energy, which indicates that the processes of their formation are spontaneous, and hence thermodynamically favorable. Meanwhile, the thermodynamic stability of the formed Al13monosaccharide complex highly depends upon the active centers of the monosaccharides and the Al13 2+ complex, which participate in the chemical coordination process, and the nature of the monosaccharide. The following complexes are identified as the most thermodynamically stable, since they exhibit the highest negative values (with regard to every considered monosaccharide as the   Figures S3 and S4. According to the results ( Figure S3), the standard entropy change is negative for all the processes, which can be explained by the fact that the spontaneous association of such reactants as monosaccharide and the Al 13 2+ complex gives the rise to a more compact/organized structure [82].
Concerning the standard enthalpy change of the interaction between the monosaccharides and the Al 13 2+ complex, it is highly negative with regard to every investigated process, which results in the formation of Al 13 -monosaccharide complex ( Figure S4). According to Equation (4), the most favorable condition for the formation of any compound is established if ∆H < 0 and ∆S > 0, which indicates that this process of formation can occur spontaneously at any given temperature. As shown above, the ∆S o of the processes of formation of all of Al 13 -monosaccharide complexes is negative, and the processes of their formation are exothermic (∆H o < 0). In this case, there is a competition between the entropy (the level of the disorder of the system) and enthalpy factors: The first parameter facilitates the reversible process (complex decomposition), while the latter one favors the forward reaction (complex formation) [83]. High negative values of ∆G o for all investigated processes result from the dominance of the enthalpy factor at the standard temperature, meaning that this temperature (T = 298.15 K) is low enough to facilitate the formation of the Al 13 -monosaccharide complex (T < ∆H o /∆S o ). However, at very high temperatures (T > ∆H o /∆S o ), the process of complex formation will not occur spontaneously.

Mechanisms of the Interactions
Based on the acquired results, the optimized geometric structures of the most thermodynamically stable Al 13 -monosaccharide complexes (with respect to every single monosaccharide considered as the reactant), and the determined mechanisms of their formation during flux enhancement in BF-MBR, applying prepolymerized aluminum chloride with the medium basicity, are demonstrated in Figure 4.
According to the determined mechanisms (Figure 4),  Figure 4g). The other represented complexes did not contain the newly formed intramolecular hydrogen bonds; thus, this factor was excluded as the major contributor to the stabilization of the investigated complexes. However, the role of the stabilizing factors in the formation of Al 13 -monosaccharide complexes should be the subject of a further investigation.

Statistical Analysis
Partial least squares analysis (PLS) of the thermodynamic parameters acquired during all simulations was conducted according to the following variables and response function ( Table 2). Table 2. Model inputs.

Predictors (X)
Response (Y) After the weighting at the first stage of the analysis to reduce the inherent differences of two predictor variables having the different ranges, the obtained model underwent random cross-validation in PLS, during which the dataset was divided into 20 segments, with 12-13 samples in each segment. The number of PLS factors was chosen according to the explained variance.
The output from the PLS analysis with regard to the scores and loadings plots is demonstrated in Figure 5.
According to Figure 5, the first two factors in total describe 100% of the variance in the dataset for x and y.
The scores plot demonstrates no clear clusters of similar samples ( Figure S5), which can be explained by a large number of category variables. Meanwhile, the samples, marked with the dark red oval, which are located leftmost, in the third quadrant of the scores plot (Figure 5a Consequently, the complexes formed by α-D-glucopyranuronic acid and the Al atoms no. 69, 73; 62, 37; and 69, 62 of the Al 13 2+ complex make up the majority of the most thermodynamically stable samples. According to the PLS correlation loadings plot (Figure 5e), Factor-1 clearly accounts for the standard enthalpy change (∆H o ) and the Gibbs energy change (∆G o ), while Factor-2 mainly describes the standard entropy change (∆S o ). All of the variables explain more than 50% of the variance, and hence have high importance in relation to Factor-1 and Factor-2. ∆H o has a high positive correlation with ∆G o and both are negatively correlated with ∆S o along Factor-1.
The explained variance plot (Figure 5f) indicates that the optimum number of factors is one, which provides 98.9% of the explained Y-variance. The increase in the number of factors does not significantly influence the explained variance, which was expected due to a relatively small number of variables (i.e., two variables result in two factors at most).
The analysis of the validation plot (Figure 5g) indicates that the developed model has a linear trend (R-squared = 0.99) with a good fit for the majority of the data (slope = 0.987). The model is reliable, since the value of R-squared (Pearson) is close to the R-squared correlation (0.989 vs. 0.995). The Root Mean Square Error of Cross Validation (RMSECV) and the standard error of cross-validation (SECV) are equal to 16.94 and 16.97 kJ/mol, respectively. The obtained bias of the model is insignificant (0.21), which indicates its low tendency to any over-or underestimation of the validation values.
The stability and applicability limit of the partial least squares regression (PLSR) model was checked through the case scenarios based on the addition of the proportional and additive noise, i.e., the noise which typically affects instrument amplification and the measurement signal, respectively, to both the predictors and response function, according to Equation (5) [84]: where PN is the level of proportional noise, %; AN is the level of additive noise, AN = P% 100 ·M(i, k), where P% is the level of approximate additive noise in percent; N(m,s) is the randomly distributed value, where m-mean standard deviation and s-standard deviation; M(i,k)-the real raw value, acquired during the quantum chemical and thermodynamic calculations; M new (i,k)-the value with the added tested noise.
According to the performed series of calculations, the derived partial least squares regression (PLSR) model is applicable (R-squared = 0.7) for the maximum proportional noise equal to 6% and the maximum standard deviation for the additive noise equal to 14, applied for both of the predictors (∆H o , ∆S o ) and the response function (∆G o ).
Due to the scattered nature of the scores, i.e., the lack of indication of distinctive clusters, it was decided to conduct the cluster analysis of the scores to understand if any grouping of the samples can be performed based on their similarities with regard to the specified variables.
The Y-scores for the latent variables Factor-1 and Factor-2 acquired during the PLSR analysis were used as the input data for the cluster analysis. This was done in order to screen out the noise of the raw data and use the prevailing differences of the most and least thermodynamically stable complexes. The classes were generated from the scores by applying a hierarchical complete-linkage clustering method and the squared Euclidean distance as the dissimilarity function.
The obtained results are represented in Figure 6 and Figure S6. The acquired dendrogram visualizes the clustering (Figure 6), which corresponds to the Y-scores in the four quadrants of Figure 5a-d. Six classes can be identified at a higher resolution (a relative distance of around 0.8) (Figure 6), which correspond to the six groups of the standard Gibbs energy change in the scores plot, generated during the data-driven sample grouping (Figure 5d).  Meanwhile, at a lower resolution (a relative distance of around 4.9), merely two classes can be identified ( Figure 6). The red cluster is characterized by the samples with high values of Factor-1 and generally low values of Factor-2. The corresponding samples were marked with the black circles in  Figure 5a-d). Hence, the corresponding complexes are more thermodynamically stable. The location of the scores of the certain sequence numbers, which constitute each cluster shown in the dendrogram, can be found in the scores plot with no sample grouping (Figure 5a).
The described clusters, generated by cluster analysis, represented the Al-monosaccharide complexes, which are more or less liable to decompose during coagulation and can be assigned to the categories "Less stable" and "Stable", respectively.

The Working Hypothesis, Implication of the Findings, Challenges, and Perspectives
Systematic research, conducted in the present work, develops the concept of the targeting chemical influence of the Al-based prepolymerized coagulant on the fouling propensity of mixed liquor in the MBR and BF-MBR systems. Particularly, this study focuses on the investigation of the thermodynamic patterns and the mechanisms of the removal of the carbohydrates, which were proven to be the major foulants, deteriorating the filtration performance of MBR and BF-MBR.
The need for this approach was emphasized by the latest work of this research team [31], which, as well as most of the other previous studies, dealt with the investigation of mechanisms of fouling mitigation via coagulation in MBR and BF-MBR in relation to the substrate, with no focus on the chemical fouling mitigation based on the carbohydrate composition [85][86][87].
Based on the analysis of zeta potential profiles and the intrinsic coagulant charges, Kulesha et al. [31] found out that the principal mechanism of chemical flux enhancement applying inorganic coagulants in BF-MBR is adsorption/charge neutralization. According to Stumm and Morgan [88], the charge neutralizing process accompanies the fixation of the multivalent cations onto the deprotonated group of macromolecular sol (carbohydrates, proteins and their derivatives). This fixation is an electrostatic or a chemical interaction. Meanwhile, the chemical coordination is of considerable significance in coagulation and flocculation reactions between multivalent cations and the macromolecular sols, producing soluble and insoluble complexes, which may have various extents of cross-linkage [88].
The underlying hypothesis of the present research is that purely chemical factors must be considered, apart from the electrostatic double-layer interactions, in order to comprehensively explain the coagulation processes occurring during the addition of the chemical flux enhancers in BF-MBR [88]. No studies on MBR/BF-MBR, which proved/refuted this hypothesis, are found.
The current study demonstrates that the investigation of the chemical coordination processes is among the key approaches to the comprehension of the fouling mitigation mechanisms in BF-MBR. The results on the reduction of the content of soluble microbial products (SMPs) after chemical dosing, acquired in our previous work [31], as well as in the studies by the other research teams [27,85,89], can be satisfactorily interpreted with regard to chemical interactions, described in the current work, while the electrosurface phenomena, described by the theories of double layer and electrostatic interactions, are not applicable to the explanation of the interactions between the solutes/macromolecular sols and the soluble hydrolyzed coagulant species.
According to Wu and Huang [90], SMPs exhibit a higher contribution to membrane fouling in MBR than the colloidal matter. Meanwhile, carbohydrates were found to be the dominant species in SMPs [91]. The carbohydrate fraction of the SMPs highly correlates with fouling in contrast to the protein fraction.
Specifically, carbohydrates were found to be responsible for the reversible, irreversible, and irrecoverable fouling, while the proteins demonstrated occasional or no quantitative relationship with the fouling extent [34,90,[92][93][94]. On the other hand, in the studies by Choi and Ng [95] and Zhou et al. [96], the proteins exhibited a high fouling propensity in the submerged MBR. Hence, the role of the proteins in membrane fouling still remains controversial, and needs further investigation.
The general mechanisms of fouling mitigation in BF-MBR were proven to be almost the same as in the MBR systems [31], which allows for assuming that the investigated content of the mixed liquor in MBR, represented by different research groups [50,51,[53][54][55], is also valid for the BF-MBR.
Due to the lack of the reference data on the thermodynamic parameters of formation of the complex in comparison to Al 13 7+ Keggin complex, and proves the hypothesis that the latter complex in its usual form is less likely to participate in any reactions [46]. The results acquired during the testing of different limit noises for the applicability of the derived PLSR model indicate that this model can be used, not merely for the computational results, but also for the experimentally-obtained data with the maximum noise, which affects the instrument amplification, 6%, and the maximum standard deviation for the noise, which affects the measurement signal, equal to 14, imposed on all the investigated variables (∆H o , ∆S o , ∆G o ).
The performed PLSR and cluster analyses are the bases for developing a classifier that enables a continuous discrimination of the thermodynamically stable and less stable Al 13 -monosaccharide complexes, based on the coagulant type and the content of the monosaccharide species in mixed liquor in the separation chamber of the MBR/BF-MBR system, which will considerably contribute to plant surveillance and process control.
The concept, presented in the current study, provides a physicochemical and statistical verification of the reasons for the observations on the SMPs removal, which were obtained during the filtration tests in BF-MBR, conducted earlier [31].
The developed approach will help to reduce the input (time and resources), required to test the coagulants as the fouling reducers in MBR and BF-MBR systems. It addresses the questions "What to expect from the mixed liquor system with the chemical dosing?" and "How to reduce the mixed liquor fouling potential?" The quantum chemical and thermodynamic calculations, followed by the multivariate chemometric analysis, allow the researcher to predict the efficiency of the coagulants with regard to SMPs removal in MBR and BF-MBR, and to select the most efficient compound without conducting numerous experiments. The thermodynamic stability of the investigated Al 13 -monosaccharide complexes is directly related to their fouling potential, since if the complex is more liable to decompose during/after coagulation, more severe membrane fouling is expected due to the released SMPs [97]. Hence, the present study introduces the relationship-fouling as a function of the standard Gibbs energy change-which is able to provide the best solution for the selection of the chemical fouling reducer in MBR and BF-MBR at the lowest cost.
The results of the present work can be used for the development of a sensor for the prediction of the flux enhancement efficiency based on the monosaccharide content of mixed liquor in BF-MBR.
Further study is also foreseen in the investigation of the role of the stabilizing factors in the formation of Al 13 -monosaccharide complexes.
Besides, the thermodynamic properties of the interaction of the monosaccharide-amino acid network assemblages with the polyaluminum hydrolysis species should be investigated.

Conclusions
The developed computational thermodynamic-multivariate chemometric approach to the assessment of the chemical flux enhancement of the coagulants with regard to MBR and BF-MBR advances the field of fouling mitigation, and provides the thermodynamic and statistical tools for the understanding of the underlying mechanisms and the prediction of the fouling mitigation efficiency.
The case study, presented in the current work, focuses on the flux enhancement efficiency of the Al 13 2+ complex, which is the main hydrolysis species of the medium basicity prepolymerized aluminum chloride, proven to be the most efficient coagulant with respect to fouling mitigation in BF-MBR. The mechanisms of the formation of the most thermodynamically stable Al 13 -monosaccharide complexes were defined based on the results of the quantum chemical and thermodynamic analyses. The thermodynamic stability of the formed Al 13 -monosaccharide complex was found to be highly dependent upon the active centers of the monosaccharides and the Al 13 2+ complex, which participate in the chemical coordination process, and the nature of the monosaccharide. The derived PLSR model can be used for computationally and experimentally acquired thermodynamic parameters, being applicable (R-squared = 0.7) for the maximum proportional noise equal to 6%, and the maximum standard deviation for the additive noise equal to 14, applied for both the predictors (∆H o , ∆S o ) and the response function (∆G o ).
The results of the performed PLSR and cluster analyses are the basis for developing a classification model for the continuous discrimination of the thermodynamically stable and less stable Al 13 -monosaccharide complexes, based on the coagulant type and content of the monosaccharide species in mixed liquor, which will substantially contribute to the automated process control of the MBR/BF-MBR systems.