We employ statistical learning to present a principled framework for the establishment of quantitative structure-property relationships (QSPR). We focus on property predictions of industrial polymers formed by multiple reagents and at varying molecular weights. We develop a theoretical description of QSPR as well as a rigorous mathematical method for the assimilation of experimental data. Results show that our methods can perform exceptionally well at establishing QSPR for glass transition temperature and intrinsic viscosity of polyesters. Keywords: QSPR, polymer, polyester, glass transition temperature, intrinsic viscosity, statistical learning 1 1 Introduction Material design is an important aspect of modern Chemistry and Chemical Engineering, especially for industrial applications [1]. Material design is particularly challenging when it comes to the design of polymers [42, 2]. Due to the complexity of modern polymer materials, the variety of available resources, and a high sensitivity on the raw materials, the design pipeline of polymers often starts with computational studies that can narrow down a small set of candidate compositions which, after appropriate processing, may give rise to materials of the desired properties [29]. Due to time and cost constraints, such studies most often rely on quantitative structure-property relations (QSPR) that are either well-established or tentative [43]. By QSPR we mean a mathematical relationship between material properties of interest and a set of structural descriptors of the molecules of the target materials desired to have such properties [32, 13]. QSPR methods are flexible as to what descriptors and mathematical models are used, but in general, QSPR maintain the formulation of their relationships as a function that maps a set of molecular descriptors to one or more numerical properties of interest [10, 31]. To date, the process of establishing a QSPR is largely achieved empirically [38]. In particular, to establish a QSPR, first a suitable set of data relevant to the targeted properties and materials needs to be identified based on expert opinion and manual curation of historical measurements. Then, an appropriate set of mathematical relations that summarizes the data needs to be identified. Subsequently, values of various unknown parameters in these relations are adjusted during parameter fitting. Finally, the pipeline is terminated with a validation stage in which the developed QSPR is tested or cross-validated against real-life data outside these used in its development. In practice, several mathematical methods are employed in each of the stages and these involve a combination of heuristics derived from experience [14] ad hoc algorithms [37] or, more recently, Artificial Intelligence implemented via Machine Learning [16]. Like most data-driven methods [18], the performance and practicality of the resulting QSPR depends critically not only on the quality of the supplied data, but also on the mathematical methods employed in its establishment [26, 39]. For this reason, it is important that all stages of QSPR development are carried out within the same framework and under a single computational process. For this, a unifying framework of data representation, structure-property relation, and parameter estimation via assimilation of experimental measurements is required. Such a framework needs (i) to be flexible enough to support arbitrarily complex compositions (for instance materials created by more than a few reagents) and (ii) to allow for parameter estimation that is computationally efficient and robust even in cases of an excessively high number of involved parameters. Here, we present such a unifying framework for QSPR. Our framework builds on chemical intuition to motivate mathematical forms for the QSPR functions and applies the principles of Statistical Learning [44, 30] to adjust the values of the unknown parameters in them. As a result, we present a comprehensive mathematical framework that can speed up QSPR development and increase its robustness. Built on statistical principles, our framework does not only provide property predictions but also estimates the uncertainty associated with them, making it particularly useful in cases of limited, corrupted, or highly complex data. The remaining of this study is structured as following. In Methods we first set basic terminology, describe our data representation scheme, QSPR forms, and parameter estimation methods. In Results we demonstrate the use of our framework by key applications in polyester materials. Finally, in Discussion we present an overview of the advantages and disadvantages of our QSPR framework. 2 Acids Reactants Abbreviations Types Functionalities Functional groups Cap acid - V 1 f1 Terephthalic acid TPA A 2 f2, f3 Isophthalic acid IPA A 2 f4, f5 1,4-Cyclohexanedicarboxylic acid 1-4-CHDA A 2 f6, f7 Adipic acid ADP A 2 f8, f9 Sebacic acid SA A 2 f10, f11 Hexahydrophthallic anhydride HHPA A 2 f12, f13 Trimellitic anhydride TMA A 3 f14, f15, f16 1,12-Dodecanedioic acid DDDA A 2 f17, f18 Azelaic acid AZL A 2 f19, f20 Alcohols Reactants Abbreviations Types Functionalities Functional groups Cap alcohol - V 1 f1 1,6-Hexane diol HDO A 2 f2, f3 Neopentyl glycol NPG A 2 f4, f5 Trimethylolpropane TMP A 3 f6, f7, f8 1,4-Cyclohexanedimethanol 1,4-CHDM A 2 f9, f10 1,3-Cyclohexanedimethanol 1,3-CHDM A 2 f11, f12 2,2,4,4-tetramethyl-1,3-cyclobutanediol TMCD A 2 f13, f14 2-Methyl-2,4 Pentanediol MPD A 2 f15, f16 4,8-Bis(hydroxymethyl)tricyclo[5.2.1.02,6]decane TCDDM A 2 f17, f18 Table 1: List of the acid and alcohol reactants we consider in our polyesters. Types A and V distinguish our reactants into actual and virtual ones. The former are those used in the preparation of a polyester; while the latter are those introduced for computational purposes only. 2 Methods 2.1 QSPR schemes Each polyester resin created has its own numerical properties, such as glass transition temperature or intrinsic viscosity, and its own composition of raw ingredients such as moles of reactants. Here, our goal is to develop a QSPR scheme that uses polyester composition to predict polyester resin properties. 2.1.1 Polyester ingredients In this study, we consider QSPR developed from polyester data that consist of compositions and experimentally measured properties. For concreteness, we consider only polyesters formed by the acid and alcohol reactants shown in table 1. These are polyfunctional monomers that react to form acid-alcohol dimers which, in turn, are our polyester’s repeat units. We classify our reactants as actual and virtual. The former are those used in a polyester’s preparation; while, the latter are those needed in our framework only for computational reasons. In particular, we include acid and alcohol caps which are monofunctional virtual monomers that can simulate the formulation of finite molecules and so to allow for the study of polyesters at varying molecular weight including those at the low end [14]. To represent each functional group of our reactants, we use the scheme shown in table 1. This scheme uses fi with subscripts i= 1, . . . , I to represents acid groups and fjwith superscripts j= 1, . . . , J to represent alcohol 3 groups. In this particular study I= 20 andJ= 18 . Similarly, we use h1:Iandh1:Jto denote the moles of each group present in a polyester. For the actual monomers, this is obtained by the moles of the corresponding reactant in a polyester multiplied by its functionality. To obtain the moles of each cap reactant, we assume that the moles of every acid and every alcohol group in a resin are equal, so a polyester’s formation terminates when both its acid and alcohol free groups are exhausted [14]. Our assumption results in the formulas h1= max 0,X j̸=1hj−X i̸=1hi , h1= max 0,X i̸=1hi−X j̸=1hj . Further, we use fj ito denote a bond formed in a polyester by the reaction of groups fiandfj. Similarly, we use hj ito denote the bond’s moles within a polyester. Assuming no direct bond formation between cap monomers, we readily obtain h1 1= 0. Further, we assume well-mixed conditions and so no preferential reactions between the groups forming each actual reactant. This results in moles of bonds obtained by hj i=hihj P ihi=hihj P jhj. Here, the denominators are equal to each other due to the inclusion of caps as described above. Because our polyesters are formed by polycondensation reactions, we use the formula mj i=mi+mj−mcondensate to obtain a weight mj icharacteristic of the bond fj i. For groups of actual reactants, miandmjare the molecular weights of the monomers carrying the groups fiandfjdivided by their functionalities. For groups of virtual reactants, these are by definition m1=m1=mcondensate where mcondensate is the weight of the molecule lost from the bond reaction. For our polyesters this is a single water molecule. Next, using the bond weights mj iand a polyester’s mole composition hj i, we calculate the weight fraction ωj i of our bond ingredients by ωj i=mj ihj iP iP jmj ihj i. Our QSPR scheme, detailed below, models the properties of a polyester material as a function of characteristic contributions by its ingredients. To represent contributions, we denote with pj ithe contribution to a polyester material’s property caused by the bond fj i. We build the contributions of each ingredient pj iusing characteristic primitives qj ias follows pj i=pref·exp qj i . (1) Here, prefis a reference value that allows for matching of the units and the scaling factor exp qj i to adjust for ingredient-specific characteristics. To encode structural descriptors in the selection of p1:J 1:Ior their primitives q1:J 1:I, we employ the simplified molecular-input line-entry system (SMILES) which is a standard way to turn a molecular representation into processable strings for chemoinformatics analysis [40]. Specifically, we denote with σj ithe SMILES string that is characteristic of the bond fj i. We obtain our SMILES σj iusing the string corresponding to the resulting dimer 4 created by the specified bond-forming reaction. We note that for an actual reactant monomer combined with a virtual monomer, such for in fj 1orf1 i, we ignore the virtual monomer and use simply the SMILES string corresponding to the actual monomer alone. Also we leave undefined any SMILES capturing virtual-to-virtual reactions, such as f1 1, since, as we mention earlier, we exclude such ingredients from our compositions. To generate our strings, we use the cheminformatic methods in the software package RDKit [23]. 2.1.2 Polyester properties Mathematically, our QSPR is a function that maps ingredient contributions, encoded via p1:J 1:I, and ingredient compositions, encoded via ω1:J 1:I, to resin properties, denoted with Ppredicted. We represent this relationship with a function Ppredicted=G(ω1:J 1:I, p1:J 1:I). Such function corresponds to a mixing rule [39] which we model by a non-linear average G ω1:J 1:I, p1:J 1:I =g−1 X iX jωj ig(pj i)  (2) induced by an invertible function g(p). Here, g−1(p)denotes the inverse of g(p)which satisfies g−1(g(p)) = p. According to our formulation, the contribution pj iis the property of a polyester material exclusively consisting of the ingredient fj i. 2.2 Statistical learning Given a mixing rule, induced by an appropriate function g(p), our QSPR contains uncharacterized ingredient contributions p1:J 1:I. Since, unlike ω1:J 1:I, these are not polyester specific, their values can be obtained via formal parameter estimation. From now on, we employ Statistical Learning [18, 30, 20] to describe a mathematical method that uses experimentally obtained measurements to derive values for p1:J 1:Ior, equivalently, values for their primitives q1:J 1:I. For clarity, we consider the assimilation of experimentally acquired data and use subscripts r= 1, . . . , R to label the individual polyester materials in our dataset. Our dataset contains Rpolyester assessments and consists of property measurements, denoted with Pmeasured r , and reactant moles, denoted with (h1:I)rand(h1:J)r. To denote the entire collection of measured properties we use Pmeasured 1:R . To process our data, we follow the Bayesian methodology [30, 36, 24]. This requires a likelihood that assumes values for the unknowns q1:J 1:Iand links predicted with measured properties. In addition, it requires a prior on the unknowns q1:J 1:Ithat penalizes nonphysical values. Below, we detail our choices. 2.2.1 Likelihood Our likelihood is a probability distribution that penalizes deviations of property predictions Ppredicted r from property measurements Pmeasured r based on assumed values of the primitives q1:J 1:I. Because, in this study, we consider only properties such as glass transition temperature or intrinsic viscosity, which are naturally positive quantities, we employ a likelihood in the form Pmeasured r |τ, q1:J 1:I∼LogNormal Ppredicted r , τ . (3) Here, Ppredicted r is the property prediction resulting from the primitives q1:J 1:Ivia eq. (1) and the mixing rule in eq. (2). Our probability distributions, such as LogNormal , are parameterized as detailed in appendix B. 5 2.2.2 Priors Our priors are probability distributions that penalize nonphysical values for the unknown parameters in the likeli- hood. For this purpose, on the likelihood’s precision τwe use τ∼Gamma ( ϕ, τref/ϕ), which excludes negative values. Our choice is conditionally conjugate to eq. (3), which allows for efficient computational implementation of our method as described below. Finally, on the primitives q1:J 1:Iwe apply a joint multivariate normal prior [8] q1:J 1:I∼NormalJ I(0, S) and use the covariance matrix Sto incorporate structural information represented by the SMILES σ1:J 1:I. Specifically, we use chemical similarity [4], encoded via the ingredients’ SMILES, to setup a matrix Ssuch that it balances an ingredient’s contribution relative to another’s. Our covariance matrix Sis a square positive definite matrix with entries Sj,j′ i,i′that are defined via the sum of a white noise and a square exponential kernel [6, 18, 30, 41]. This results to Sj,j′ i,i′=λ2 a δj,j′ i,i′+ (1−a) exp −1 2 1− T M σj i ,M σj′ i′ ℓ 2   (4) where a, ℓ,andλare tunable parameters and δj,j′ i,i′is the Kronecker delta function. In turn, our covariance kernel is activated by the Tanimoto coefficient [3] and the Morgan fingerprint [7] denoted with T(µ, µ′)andM(σ), respectively. We compute these via RDKit [23]. Our choice of the covariance in eq. (4) preserves the comparison coefficient of similarity or difference but can scale the relationship to reflect how the similarity in structures relates to the similarity in measurements [41, 11]. In essence, as we explain below, with our combination of Morgan fingerprint and Tanimoto coefficient, we have a comparison metric that quantifies the intensity of similarity or dissimilarity between two SMILES strings [3]. We use this in eq. (4) as a basis to make inferences on the correlation of the contributions associated with our QSPR ingredients. Put simply, our covariance Sfits similar contributions pj ito ingredients fj iwith similar SMILES σj i and vice versa dissimilar contributions to ingredients with dissimilar SMILES. Chemical information on the Morgan fingerprint µ=M(σ)comes in the form of a bit vector µ= (µ1, µ2, . . .) with each bit µkindicating whether an attribute is present or not in that spot. To calculate the Morgan fingerprint we also need to specify the radius (number of atoms used to create a mapping) and the nBits (the number of unique combination pathways to represent features). Once a radius and a bit size are chosen, the Morgan fingerprint is created by (i) identifying each atom in the assembly, (ii) updating each atom’s identifiers based on its neighbors (to the size of the radius), (iii) removing duplicates, and (iv) collapsing down into a single-bit vector. With the Morgan fingerprint representation M(σ)of a SMILES σin hand, we then define a metric for comparison of our SMILES. We opt for the conventional choice of the Tanimoto similarity coefficient [34]. This is defined by T(µ, µ′) =µ·µ′ ||µ||2+||µ′||2−µ·µ′ forµ=M(σ)andµ′=M(σ′), note that for self-comparison we get T(µ, µ) = 1 . 6 (ω1:I1:J)r p1:I1:J prefq1:I1:Jλ ℓ σij i=1,...,I j=1,...,Jresin data ingredient data primitivescontributionsS1:I,1:I1:J,1:J ττref φ r=1,...,RPrpredicted PrmeasuredαFigure 1: Schematic representation of our learning framework: Our QSPR modeling scheme includes ingredient primi- tives q1:J 1:I, ingredient contributions p1:J 1:I, weight fractions (ω1:J 1:I)1:R, and property measurements Pmeasured 1:R . The latter two are polyester specific; while, the former two are not. The SMILES of our ingredients σ1:J 1:Iare used to generate a covariance matrix S1:J,1:J 1:I,1:I, with hyperparameters λ, α, ℓ that controls the similarity between the fitted primitives q1:J 1:I. In turn, fitted primitives q1:J 1:Idetermine the contributions p1:J 1:Iand, via the mixing rule, eventually the predicted polyester properties Ppredicted 1:R . Following the common convention: we denote random variables with circles, deterministic variables with diamonds, and we leave open hyperparameters with set values. Measurements are shaded and arrows indicate the dependencies between our variables. Finally, we highlight in color the dependencies mediating our priors. 2.3 Computational inference Our statistical learning framework is summarized in fig. 1 and it is fully detailed in appendices A.1 and A.2. We use our scheme to conclude on unknown parameter values via its marginal posterior probability distribution P q1:J 1:I|Pmeasured 1:R . Using Bayes’ theorem [36], the posterior is formally derived by P q1:J 1:I|Pmeasured 1:R ∝ P q1:J 1:IZ τdτP(τ)Y rP Pmeasured r |τ, q1:J 1:I , (5) where P(q1:J 1:I)andP(τ)are our prior distributions and P Pmeasured r |τ, q1:J 1:I our likelihood. Because of the complexity of eq. (5), we cannot obtain a closed form. Instead, we characterize our posterior using Markov chain Monte Carlo (MCMC) sampling by drawing samples of q1:J 1:I|Pmeasured 1:R [30, 35, 25]. Because our model attains the form of a latent Gaussian model, that is a multivariate Gaussian prior and a complicated likelihood function that ties the unknown variables to the observed data, we develop our MCMC scheme by means of an elliptic slice sampler [27]. 2.4 Data acquisition To generate the results shown later in this study, we relied on polyester data. Our data are acquired by poly- condensation and characterization procedures as described next. 2.4.1 Reagents Isophthalic acid (IPA), terephthalic acid (TPA), 1,4-cyclohexanedicarboxlyic acid (1,4-CHDA), 1,4- cyclohexane dimethanol (1,4-CHDM), 1,3-cyclohexane dimethanol (1,3-CHDM), 2,2,4,4- tetramethylcyclobutanediol (TMCD), hexahydrophthallic anhydride (HHPA), and neopentyl glycol (NPG) were provided by Eastman Chemical Company. 7 Adipic acid (AD), sebacic acid (SA), trimellitic anhydride (TMA), 1,12-dodecanedioic acid (DDDA), azelaic acid (AZL), 2-methyl-1,3-propanediol (MP Diol), trimethylolpropane (TMP), 1,6-hexane diol (HDO), and 4,8- bis(hydroxymethyl)tricyclo[5.2.1.02,6]decane (TCDDM) were purchased from Sigma Aldrich. Aromatic 150 was purchased from ExxonMobile and Fascat 4100 was purchased from PMC Organometallix Inc. 2.4.2 Polyester resin synthesis The polyols were produced using either a resin kettle reactor setup via solvent-assisted polycondensation or a resin rig reactor setup via melt polycondensation, both of which were controlled with automated control software. The solvent-assisted resins were produced on a 3.5 mole scale using a 2 L kettle with overhead stirring and a partial condenser topped with total condenser and Dean Stark trap. Approximately 10 wt % (based on reaction yield) azeotroping solvent of high boiling point (A150 and A150ND) was used to both encourage egress of the water condensate out of the reaction mixture and keep the reaction mixture viscosity at a reasonable level using the standard paddle stirrer. Chemical reagents were added to the kettle, which was then completely assembled. The Fascat 4100 (monobutyltin oxide) catalyst was added via the sampling port after the reactor had been assembled and blanketed with nitrogen for the reaction. Additional A150/A150ND solvent was added to the Dean Stark trap to maintain the ∼10 wt % solvent level in the reaction kettle. The reaction mixture was heated without stirring from room temperature to 150 °C using a set output controlled through the automation system. Once the reaction mixture was fluid enough, the stirring was started to encourage even heating of the mixture. At 150 °C, the control of heating was switched to automated control and the temperature was ramped to 230 °C over the course of 4 h. The reaction was held at 230 °C for 1 h and then heated to 240 °C over the course of 1 h. The reaction was then held at 240 °C and sampled every 1-2 h upon clearing until the desired acid value was reached. The∼90%-solids resins were ground into 6 mm pellets and thoroughly dried in a vacuum oven at 150 °C for 24 h prior to characterization. The melt polycondensation resins were produced on a 0.5 mole scale. A 500 mL, one-neck, round- bottom flask was carefully charged with all chemical reagents and Fascat 4100 (monobutyltin oxide) catalyst. The flask was equipped with a polymer head adapter with stainless steel mechanical stirrer and securely clamped to the polymerization rig. To the polymer head, a distillation side arm and Erlenmeyer flask were attached. The automated-controlled vacuum system was attached to the flask side arm to allow for a reduction in pressure of the reaction vessel. The Belmont metal bath was preheated to 20 °C above the recipe starting temperature (180 °C). The apparatus was subjected to two iterations of a nitrogen ( N2) purge to remove oxygen and then dunked into the metal bath to begin. The flask was held at 180 °C for 10 minutes to melt the starting materials and then stirring was started to encourage even heating on the mixture. The flask was heated to 240 °C over 4 hours and then held there for an additional hour. Pressure in the reaction flask was reduced to 1.5 torr over 45 minutes and then the reaction was held at 1.5 torr until the final acid value was reached. Dry ice was used to ensure that the solvent traps were sufficiently cold to prevent any solvent/organic matter from going to the vacuum pump. After completion, the flask was slowly brought back to atmospheric pressure and removed from the hot metal bath. Upon solidification, the polymer was pulled from the round bottom flask by partially melting the edges, then the glass flask was broken with a hammer to give the solid polymer “lollipop” on the stir rod. The polymer was cooled in dry-ice, removed from the stir-rod, and ground into 6 mm pellets prior to characterization. 8 2.4.3 Polyester resin characterization The Inherent Viscosity ( [η]) of all polymers was determined in 0.5 wt% PM 95 (60/40 phenol/1,1,2,2- tetra- chloroethane) solution at 25 °C. Molecular weights were determined by gel permeation chromatography (GPC) with 95/5 methylene chloride/HFIP mobile phase and calibration curves for polystyrene standards. Monomer composition was determined via gas chromatography (GC) hydrolysis. The glass transition temperature (Tg)was determined using differential scanning calorimeter (DSC) at 20 °C / min. ramp rate with a N2sweep. The Tg was based on second heat thermograms. 3 Results In this section, we demonstrate how our statistical learning framework applies in the establishment of QSPR schemes for polyester resins. Our raw data assess polyesters created by 9acid and 8alcohol reactants with different functionalities resulting in I= 20 andJ= 18 different functional groups, forming 360different bonds, as shown in table 1. The calculated Tanimoto coefficients and associated prior covariance matrix Sis detailed in appendix D and the sensitivity of our method to the hyperparameters is detailed in appendix E. For each polyester resin, our data contain reactant mole amounts h1:I, h1:Jand measurements of glass transi- tion temperature Tg, intrinsic viscosity [η], and number average molecular weight Mn. From these, Tmeasured g and [η]measuredare our targeted properties. Overall, we use R= 137 resin data points with Tmeasured g measurements, R= 79 resin data points with [η]measuredmeasurements, and R= 79 resin data points with both Tmeasured g and [η]measuredmeasurements. For all cases shown below, we obtain the maximum a posteriori estimator ˆq1:J 1:Iof our primitives by selecting the MCMC sample with the highest P(q1:J 1:I|Pmeasured 1:R )computed via eq. (5) [30, 6, 36, 18]. This estimator represents the single best choice of parameter values generated by our learning framework. Using this estimator, we obtain our contributions ˆp1:J 1:Ivia eq. (1) and our QSPR predictions ˆPpredicted r for the rthresin using our mixing rule ˆPpredicted r =G (ω1:J 1:I)r,ˆp1:J 1:I . 3.1 Glass transition temperature First, we demonstrate how our method applies in the case of establishing a QSPR for Tg. For our QSPR development, we use pref= 65.30◦C and g(p) = 1 /p. Specifically, fig. 2 shows maximum a posteriori predictions ˆTpredicted g and compares them with Tmeasured g . The credible intervals result from taking the whole set of sampler predictions calculated under our qI:J 1:Ivalues and summarizing them by their 5thand 95thpercentiles. For comparing our (ˆTpredicted g )1:Rwith (Tmeasured g )1:Rwe use the coefficient of determination [17]. This is denoted with ˆR2and detailed in appendix C.1. The resulting values of ˆR2quantify how well properties from the data are reproduced by our QSPR. In particular, ˆR2= 0 indicates no ability to make a prediction and ˆR2= 1 indicates a perfect replication of the data [9]. We emphasize that our learning scheme does not rely on ˆR2, and therefore ˆR2is not used in our learning process. We employ it here only to facilitate a presentation of how well our model has been able to match the data it was supplied with. To further assess the performance of our framework, we explore leave-one-out cross-validation of our model by holding one of the resin data points for our Tgdata out of our learning process and then using the contribution pre- dictions of our model to generate predictions for the Tmeasured g of the withheld resin. Here we report the experimen- tal data values as well as the 5th, 50th, and 95thpercentiles denoted ((Tpredicted g )5%,(Tpredicted g )50%,(Tpredicted g )95%) 9 Figure 2: Quantitative comparison of our QSPR method using the maximum a posteriori estimator ˆTpredicted g , plotted against the experimentally determined Tmeasured g from the resin data. In this figure, points represent predicted and measured values. Points along the dashed line represent the perfect agreement between our prediction and measurements data and any deviation from this line is factored in the ˆR2value. For this figure, we get the value ˆR2=.96. Figure 3: Leave-one-out cross-validation comparison of Tpredicted g withTmeasured g : Histograms summarize suggested prop- erty values and vertical lines compare with their measured ones. For a numerical comparison see table 2 below. as well as the error of our percentiles compared to the measured value Terr5%g, Terr50%g , Terr95%g see appendix C.1 for the formulas used to calculate the error. As the selected results in fig. 3 indicate, the experimental values typ- ically fall between the ((Tpredicted g )5%,(Tpredicted g )95%)interval and often get close to our median (Tpredicted g )50% prediction. To demonstrate how effective our approach is, we test it against a well-validated method on the same data for Tmeasured g . Specifically, we compare our methods against the formulation of Bicerano [5] as shown in fig. 4. We also note that with Bicerano’s prediction method, we found it necessary to correct the infinite molecular weight predictions with additional molecular weight data, namely the number average molecular weight Mn. Our method forTpredicted g does not require this information, but is limited to predictions for resins in the same molecular weight neighborhood. For further explanation of the methods used to generate the results, we compared our methods against see appendix F. 3.2 Intrinsic viscosity Because previous studies [33, 12] suggest a strong dependence of [η]on molecular weight that is captured by the empirical Mark–Houwink equation [39], we use our QSPR scheme to model improvements to the Mark–Houwink 10 Tmeasured g Tpredicted g Comparison Resin (Tpredicted g )5% (Tpredicted g )50% (Tpredicted g )95% Terr5%g Terr50%g Terr95%g 1 89.52 86.11 89.74 92.13 -3.80 0.25 2.92 2 90.42 89.16 90.47 92.10 -1.39 0.06 1.85 3 55.58 53.24 56.60 59.72 -4.20 1.84 7.45 4 61.80 59.11 63.46 67.47 -4.35 2.68 9.17 5 87.70 84.02 87.02 91.30 -4.19 -0.77 4.11 ◦C◦C◦C◦C % % % Table 2: Leave-one-out cross-validation comparison of Tpredicted g withTmeasured g results: A comparison of our predictions with experimental measurements reveals good agreement. This table shows experimental values against characteristic quantiles values of our predictive probability distribution. Characteristically relative error sticks below 10 percent. Figure 4: Quantitative comparison of QSPR developed using our framework via the maximum a posteriori estimator ˆTpredicted g and the Tpredicted g produced by Bicerano’s methods against the experimentally determined Tmeasured g . In this figure, points represent predicted and measured values. Points along the dashed line represent the perfect agreement between our prediction and measurements and any deviation from this line is factored in the ˆR2, R2values. For this figure, we get the value ˆR2=.96and for Bicerano’s method R2=.83. 11 Figure 5: Quantitative comparison of our QSPR method using the maximum a posteriori estimator [ˆη]predicted, plotted against the experimentally determined [η]measuredfrom the resin data. In this figure, we see blue points that represent the value for each resin data point r= 1, . . . , R from both the data and our prediction. Points along the orange dashed line represent the perfect agreement between our prediction and the experimental data and any deviation from this line is factored in the ˆR2value. For this figure, we get the value ˆR2=.71. equation instead of modeling [η]directly. Specifically, to develop a QSPR for [η]we target an artificial numerical property given by the ratio P=[η] κMαn that is insensitive to molecular weight. Subsequently, we convert our predictions PpredictedforPto predictions [η]predictedfor[η]according to [η]predicted=κMα nPpredicted. Here, κMα nis the Mark–Houwink prediction [39] with the parameter values κ, αadjusted to best match our data set. For our QSPR development, we use pref=.0197 (dL/g) and obtain the mixing rule also by g(p) = 1 /p. As earlier, for comparing our [ˆη]predicted 1:R with[η]measured 1:R we use the coefficient of determination [17], denoted with ˆR2and detailed in appendix C.2. We compare the prediction of our QSPR method for [η]predictedagainst the measurements. As earlier, we show the maximum a posteriori predictions [ˆη]predictedand compare them with [η]measuredin fig. 5. We explore leave-one-out cross-validation of our model by holding one of the resin data points for our [η] data out of our sampling and then using the repeat unit property predictions of our model to generate a col- lection of [η]predictedvalues of the withheld resin. As earlier, here we report the 5th, 50th, and 95thpercentiles denoted ([η]predicted 5%,[η]predicted 50%,[η]predicted 95%)as well as the error of our percentiles compared to the measured value ([η]err 5%,[η]err 50%[η]err 95%)see appendix C.2 for the formulas used to calculate the error. 3.3 Multiproperty QSPR Having demonstrated how our method performs in the simple context of a single property, either Tgor[η], we now move to the more demanding setting, where with our QSPR we target both properties simultaneously. Our 12 Figure 6: Leave-one-out cross-validation comparison of [η]predictedwith [η]measured: Histograms summarize suggested property values and vertical lines compare with their measured ones. For a numerical comparison see table 3 below. [η]measured[η]predictedComparison Resin [η]predicted 5%[η]predicted 50%[η]predicted 95%[η]err 5%[η]err 50%[η]err 95% 1 0.279 0.267 0.289 0.309 -4.16 3.51 10.66 2 0.321 0.318 0.338 0.361 -0.82 5.37 12.48 3 0.305 0.280 0.310 0.340 -8.21 1.54 11.28 4 0.199 0.209 0.218 0.228 5.13 9.65 14.48 5 0.479 0.391 0.428 0.466 -18.31 -10.66 -2.73 dL/g dL/g dL/g dL/g % % % Table 3: Leave-one-out cross-validation comparison of [η]predictedwith[η]measuredresults: This table shows experimental values against characteristic percentile values of our predictive probability distribution. We have a higher relative error compared with our Tgfits which is to be expected. As the selected values suggest we tend to see a good agreement of our median prediction when the value is around the average of our data but our prediction struggles as we move higher or lower. learning scheme extends naturally to this setting as detailed in appendix A.2. Namely, we now consider a bivariate property with Pr= (Tg,[η])r. Now, in fig. 7, we demonstrate how the fits work across Tgand[η]simultaneously. 4 Discussion In this study, we presented a statistical learning approach to the development of QSPR for polyester properties and, to demonstrate characteristic predictions, we assimilated data from experimentally characterized poly-condensation reactions. To model how the structural contributions of the repeat unit ingredients combined to give rise to the recorded properties we developed a novel group contribution model in the style of Van Krevelen’s foundational work [39] as well as taking a statistical learning approach to adjust the arising parameters with unknown values [30, 18]. In this study we focused on polyester resins, however our methodology is trans-formative and can be extended to other types of polyester products and even other chemical materials. Our key contribution is a method that is grounded in the statistical learning methodology as well as the recent progress of empirical methods in the field of polymers. In particular, we introduce a probabilistic formulation for the known and unknown variables involved in the development of a QSPR. This allows us to take advantage of Bayesian concepts [36, 30] and encode our entire scheme in a single posterior P(QSPR |data) . We also took the new approach of encoding the structural difference based on the molecular footprints into the prior’s co-variances of our repeat unit property predictions in the statistical learning scheme we applied to the QSPR model. From these contributions, we were able to assembly a modeling procedure that generates estimates of contributions from our polymer ingredients to give rise to accurate property prediction. We used two key properties from the experimental data and showed good agreement on both properties when fitted both apart and 13 (a)TgˆR2=.95 (b)[η]ˆR2=.67 Figure 7: Quantitative comparison of our both of our QSPR methods coupled together using the maximum a posteriori estimator ˆTpredicted g and[ˆη]predicted, plotted against the experimentally determined Tmeasured g and[η]measuredfrom the resin data. In both figures, we see blue points that represent the value for each resin data point r= 1, . . . , R from both the data and our prediction. Points along the orange dashed line represent the perfect agreement between our prediction and the experimental data and any deviation from this line is factored in the ˆR2value. For the figure, we get the values ˆR2=.95 forTg(left) and ˆR2=.67for[η](right). together. Our methods are not without their limitations, we are required to integrate new structural information into a statistical formulation. Also approach is more designer and requires more mathematical ability. To keep our method simple we ignore the cap molecules in our Tanimoto. Our method is computationally expensive compared to naive methods like a pre-trained Bicerano model, taking a couple hours on a standard research desktop with running on a single core to yield good results. And up to ten to twelve hours on the same equipment to give excellent results. However, in the context of some molecular dynamics simulations this time frame is quick. With this in mind, our main computational bottleneck is the Metropolis-Hastings scheme we use to sample our target distribution. Methods to parallelize this procedure or simplify the proposal setup could give significant speed-ups and will be explored in addition to our other extensions. In addition to computational or algorithmic improvement. Future research can try to investigate how statistical learning can be applied to non-polyester polymers. It also remains to show if statistical learning can be applied to properties beyond Tgor[η]. These may include properties like the melting point ( Tm) or equilibrium melting point ( T◦ m) [39, 19]. Acknowledgements This work is based upon raw data provided by Eastman Chemical Company, Kingsport, TN, USA. The authors thank Erik Burke for help in synthesizing a majority of the polyester resins used and acknowledge support from Eastman via grant No R011052266. Conflict of interest SM, DO, MD, and IS declare no competing interest. BA and CB declare no competing non-financial interest but competing financial interest as employees at Eastman Chemical Company. 14 Author contributions All authors conceived research and wrote the manuscript. CB contributed experimental data. DO contributed the software implementation of Bicerano’s methods. SM contributed software implementation of the developed algorithms for statistical learning and produced synthetic data. SM, DO, MD, IS developed theoretical methods. IS and MD oversaw all aspects of this work.