{'fraction_non_alphanumeric': 0.03699248713347662, 'fraction_numerical': 0.010135468223138051, 'mean_word_length': 5.020180436847103, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Recent advances in accelerated discovery through machine learning and statistical inference

Recent applications of machine learning and statistical inference provide case studies demonstrating how such approaches can accelerate the discovery process in physical chemistry and related fields. Examples discussed in this review include the introduction of automated approaches to systematically improve experimental design, increase the efficiency of computationally expensive molecular simulations, facilitate construction of predictive models for complex biological processes, and discover interparticle potentials that lead to materials which meet specified design goals. A common theme is the synergy between experiment and computation enabled by such approaches.


Introduction
Modern computation is integral to research in the physical sciences. In the present "Age of Big Data", [1][2][3] there are increasing opportunities for computational research to facilitate scientific discovery via Artificial Intelligence (AI). Historically, much of the academic work in AI was beset by disappointment, with early promises left unfulfilled, culminating in a steep decline in funding for AI research in the 1970s called the "AI winter". [4][5][6] In some ways, the overly optimistic predictions of the 1960s have still not come to pass (e.g., the sentient computer HAL 9000 from Arthur C. Clarke's Space Odyssey). [7] However, there have also been compelling successes, including the victory of Google's "AlphaGo" over a professional human opponent, that have been fueled by a resurgence of interest in AI. [8][9][10] Current work in the field is often guided by a practical data-driven goal: the ability to automatically extract specific insights from (or make predictions based on) large data sets. [5,11,12] Thus, while sentient robots may not exist, more subtle achievements have been realized including the introduction of powerful machine learning (ML) and statistical inference (SI) approaches. Such methodologies have the potential to be exploited in the physical sciences in order to 1) extract key insights, 2) make research processes more systematic, and 3) increase research efficiency. This review details recent progress towards this end.
Both ML and SI take as input an array of feature variables and attempt to fit some model to a known "training" set of data. Fitting the model involves "learning" the values of the tunable model parameters via a systematic optimization technique. [11,12] The distinction between ML and SI that we adopt in this work lies in how the model is used.
Typically, ML employs powerful black-box models to accurately predict some unknown quantity of interest that is believed to be dependent on the features. In model fitting, the training data is supplemented with prior observations of the relationship between features and the unknown quantity. For SI, the goal is typically to model the underlying process by which some set of features is generated. The process that generates the observables of interest is assumed to be representable by a probability distribution with parameters that are "learned" or "inferred".
ML and SI are not wholly new to the physical sciences. For example, widely used regression analysis methods are fundamentally ML methods. More recently though, increasingly advanced ML and SI strategies are being leveraged to solve a variety of challenging problems. In this review, we first describe the use of ML predictions to expedite the iterative experimental design process of synthesis and measurement for maximizing a quantifiable design goal. [13][14][15][16] By using ML-based methods, the results of previous experiments can be used to predict the most promising new sample to synthesize or new experiment to perform. In another study, ML is leveraged to expedite the computation of quantum mechanical (QM) forces for use in molecular dynamics (MD) simulations. [17][18][19] Prior QM calculations are used to create a predictive ML model relating configurations to QM forces.
Applications of SI in physical chemistry have also emerged to address a variety of challenges. First, we review the recent application of probabilistic statistical mechanical approaches for modeling complex biological phenomena. Application to virus models yields insight into how various residues interact, which could potentially expedite vaccine design. [20][21][22][23] Similarly, SI methods have been employed to extract intra-cellular interaction networks on the basis of DNA microarray data. [24] Finally, the review ends with discussion of two new inverse design strategies based on SI for the systematic discovery of materials which either exhibit complex self-assembly behavior [25][26][27][28][29] or possess a desired material property [30]. Both methods learn the requisite interactions via relatively simple optimization schemes.

Acceleration of conventional optimization and simulation schemes
A straightforward application of ML methodology is the prediction of a quantity of interest, F , as a function of known input conditions, X ≡ {x i }, called features; here, {...} indicates a set or vector of values of some unspecified size. Generally, F is difficult or costly to evaluate. [11,12,31] Therefore, accurate predictions are useful to save time and resources. For so-called supervised ML strategies, a predictive model, F (X|θ), is fit via flexible parameters θ ≡ {θ i } to the prior observations of the correspondence between X and F , i.e., the training data, D ≡ F (i) , X (i) . The model allows for prediction of F at X not included in the training set. Uncertainties associated with the predictions for F can be also obtained to yield a probabilistic model for the data, P (F |X, θ). The uncertainties can either be directly obtained via ML approaches (e.g., Gaussian process models) [11,12,31,32], or computed via bootstrapping [12], i.e., re-sampling the data to generate new model fits that probe the sensitivity of the model to fluctuations in D. For either approach, additional sampling both improves the prediction and reduces the uncertainty in P (F |X, θ).
Supervised ML methods essentially amount to interpolation strategies, each with different underlying assumptions regarding function form, smoothness and possibly the degree of interaction among the features. For example, the popular Gaussian process model requires a "kernel" that governs how correlations between two points in feature space decay, the choice of which influences the smoothness and uncertainty of the predicted function. [11,12,31,32] Another important consideration for ML is the construction of the features. It is oftentimes advantageous to purposefully engineer features to reflect the underlying physics of the problem, instead of requiring the model to learn these relationships via sampling. Feature engineering is non-unique and requires careful consideration of the problem at hand but can ultimately lead to better convergence of the model. [11,12] Examples of some feature choices will be discussed in the next two subsections, which are organized as follows. In Sect. 2.1, we consider how prediction of the figure of merit (or fitness) associated with a set of tunable conditions can be used to accelerate certain optimization problems. In Sect. 2.2, we discuss the use of ML techniques to approximate QM forces for use in MD simulations.

Optimization
There are myriad optimization schemes that can be used to maximize a desired property, termed the figure of merit (FOM), by tuning a set of adjustable parameters, θ. [11,12,33] Most require the use of derivatives with respect to θ. However, there are also many cases where derivatives are not easily obtained, e.g., when analytical derivatives are not available (or do not exist) and there are too many adjustable parameters for numerical differentiation to be practical. For such cases, heuristic optimization strategies such as simulated annealing or genetic algorithms have found widespread use. [33] In these methods, updates to parameters are governed by a tension between increasing the FOM and promoting a wider exploration of parameter space through stochastic processes. As a result, a large number of FOM evaluations can be required to arrive at a solution, limiting the practicality of such strategies, especially for FOM evaluations that require expensive computations or experiments. To handle such cases, ML techniques can be employed in order to generate parameter updates that are informed by prior FOM evaluations as opposed to stochastic processes that only consider possible changes to the current FOM.
Bayesian optimization (BO) represents one such approach, whereby a probability distribution, P (F |X, θ), is constructed to model the true underlying FOM function or landscape, F 0 ≡ F 0 (X), by fitting to prior observations (i.e., D above) via the tunable model parameters (θ). [31,32] P (F |X, θ) provides information on regions in parameter space that are 1) highly likely (determined via prior FOM evaluations) to possess high F values and 2) those of low certainty but large positive uncertainty in F . The former can be exploited to further refine the best-known parameters, and the latter represent profitable areas for exploration-allowing the optimization to progress to other maxima as needed. As more observations are included, progressively better predictions for F 0 (X) are generated.
After construction of P (F |X, θ), the model must be translated into a prediction for a set of parameters θ for the next optimization step. [31,32] Because this mapping is non-unique, multiple options exist for this task, each providing a particular balance of exploitation and exploration. One such choice for a recommender function, R(X|θ), is called Expected Improvement where F * is often taken as the best FOM evaluation at any prior X. [32] Only regions of uncertainty that are better than the current optimal set of parameters are included, i.e., uncertainty in the negative direction is not penalized. Including the uncertainty associated with the predictions in the recommender function has been shown to enhance the optimization efficacy. [16] As the optimization progresses, uncertainty will decrease in regions of parameter space that are sampled, perhaps directing the optimization to explore areas with greater uncertainty. BO can be generalized to suggest the most promising "batch" of samples, of size q. [31,32] In computational research, q may be the number of parallel processes that can be run, while in experiment, q might be governed by the sample auto-feeder size, for instance. One generalization adopted in the work below is called the Kriging believer algorithm, where the model is iteratively updated via predictions for F at suggested X within the same batch (i.e., "believe" the model). Once actual measurements are taken, the measured FOM values replace the predicted ones in the data set.
One recent study integrated BO into the experimental design of a NiTi-based shape memory alloy (SMA) with minimized thermal hysteresis. [13] SMAs exhibit shape memory and super-elasticity as the result of an underlying high-temperature austenite and lowtemperature martensite phase transition. The temperature at which the transition occurs depends on whether the material is being heated or cooled, yielding a hysteresis temperature gap, ∆T . Larger values for ∆T result in greater performance degradation-motivating the chosen FOM definition: F = −∆T . Optimization targeted an alloy of the form Ni 50−x−y−z Ti 50 Cu x Fe y Pd z where x, y, and z are doping elements allowed in increments of 0.1% within the constraints 50 − x − y − z ≥ 30%, x ≤ 20%, y ≤ 5% and z ≤ 20%yielding 797,504 possible combinations. The combinatorial challenge combined with the expensive FOM evaluation (experimental synthesis) called for a Bayesian approach.
BO, like all ML approaches, requires defining the features (X) to be used as input to the ML algorithm. While x, y and z are the tunable experimental parameters, their effect on ∆T is mediated by their impact on various atomic properties such as valence electron number. Therefore, the authors represented features via weighted fractions (e.g., f x = x/(x+y+z)) of the atomic properties listed in the second column of Fig. 1. Each property is a coarse-grained electronic structure measure known to influence the phase transition temperatures and hysteresis. While the ML model might eventually learn these relationships with sufficient sampling, directly encoding them is likely to improve the convergence rate of the model.
The overall iterative design strategy for alloy discovery is summarized in Fig. 1. The initial D was composed of 22 samples, which for each, differential scanning calorimetry measured ∆T . From these samples, the features were computed and used as input to the probabilistic model, from which a batch of new unexplored alloys with the highest expected improvements were drawn. The FOMs (−∆T ) associated with the new parameter combinations were experimentally determined, and these data points were then used to update the probabilistic model, from which new, more accurate, parameter combinations were extracted in an iterative fashion. After only nine iterations of q = 4 batch updates, 14 out of the 36 total new alloys posses a ∆T < 3.15K, which was the best of the initial 22 training samples. The outcome of the BO was a 42% improvement over the best initial training sample; the alloy Ti 50.0 Ni 46.7 Cu 0.8 Fe 2.3 Pd 0.2 with ∆T = 1.84K was discovered.
In conclusion, BO is a promising avenue for assisting in experimental (or computational) design. The approach is general and can be adapted to accommodate other scenarios. For example, the same authors recently leveraged BO for the accelerated search of BaTiO 3based piezoelectrics to minimize the curvature of a phase boundary. [15] Utilizing various descriptive features, prior phase diagrams, and the predicted functional form (from Landau theory) for the phase boundaries, materials with better piezoelectric properties than any in the initial training set were discovered. Beyond hard materials, the field of soft materials design could benefit from BO as 1) database knowledge is limited, and 2) experiments are often highly customized and expensive. Conveniently, many simplified theories exist for soft materials, [34,35] the predictions of which can be incorporated into the optimization as features to embed prior knowledge and improve model convergence.

Molecular Simulation
MD simulations require many interparticle force evaluations especially for systems with large numbers of particles, motivating the development of classical force fields. Such models invoke a selected functional form (e.g., a quadratic function for a chemical bond) with flexible model parameters (e.g., the associated spring constant) that are fit to reference data, which might be from either QM calculations or experiments. However, classical force fields are not always appropriate; an archetypal failure is that a harmonic spring cannot model bond-breaking processes. Ideally for such cases, QM data would be used to determine the forces for MD simulations; however, QM calculations that are sufficiently accurate to produce meaningful MD simulations are often too computationally expensive.
A ML model that is able to predict QM forces for a given particle arrangement would significantly reduce the computational expense of MD simulations with QM forces. Such a strategy would be particularly amenable to standard ML approaches due to the high degree of structural similarity between many of the configurations that arise in a given MD simulation. In general, proper sampling in MD studies results in an ensemble of configurations possessing similar structural correlations. Therefore, ML models can be employed to learn the complex relationships between QM forces and the configurations for some initial set of structures that initialize the model (i.e., the training data); the resulting model can then be used to predict QM forces. [17][18][19] The above scheme is summarized in Fig. 2. In Fig. 2a, the energy trace of a sample trajectory is shown; green regions indicate steps for which QM calculations are needed. In this case, the model is initialized using the early steps of the trajectory until sufficient sampling occurs such that an accurate predictive model can be constructed, at which point ML is used to predict the QM forces (the peach colored region). Later, a large fluctuationsignaled by the jump in the energy-produces a configuration that is outside of scope of the training data; see Fig. 2b. At this point, explicit QM calculations are performed to augment the training data until this new area is "learned" sufficiently well that the ML predictor can be used again. In the future, large fluctuations that span these two domains can be treated using the augmented ML model. Fig. 2c depicts this procedure as a flowchart. The above approach was employed to perform MD simulations of atomic Aluminum (Al) using density functional theory (DFT) to compute the QM forces for the training data. [17,18] Intuitively, the features-the known quantities from which the forces are predicted-must depend on the environment of the selected particle, in this case the other Al atoms in the simulation. In order to explicitly encode the physics of the problem, atomic coordinates were mathematically transformed into features that are invariant to any rotation, translation or permutation of identical particles. Additionally, the forces on a given Al atom should most strongly depend on its closest neighbors, with more distant Al atoms contributing less strongly. Therefore, both a damping function and a cutoff distance were included in the features to encode the primarily local effects of neighboring particles on the forces.
Inputting these engineered features into a standard ML model enabled quick and accurate force prediction. The errors in ML forces were found to be within the expected accuracy of DFT itself while reducing the computational expense by multiple orders of magnitude with respect to DFT. Such approaches seem promising and may enable QM MD simulations to reach time and length scales not previously believed to be possible. However, further research is needed to address certain open questions, e.g., how does one detect when the model is outside its regime of predictability? The approach adopted by the authors assumed any feature vector lying outside a hypercube encompassing the training set needed explicit DFT computation. While certainly true, points within the sample region could also warrant explicit computation. One possibility includes modeling not only the expected value but also uncertainty (as done in Bayesian optimization). If the uncertainty exceeds some threshold, an QM computation could be performed instead. Some recent work based on the distance between the new feature set and the closest neighbor in the training data shows promise. [19] 3 Modeling biology with statistical mechanics In Sect. 2, we reviewed use of ML techniques that prioritize the accuracy of predictions, in general leaving the relationships between the features and the predicted values obscured. In this Section, we consider applications of SI where a primary goal of the modeling is to understand the underlying physics that gives rise to a set of observables, potentially in conjunction with making accurate predictions for a quantity of interest. In the following studies, SI methods are used to construct a model, P (X), for some underlying probability distribution, P 0 (X), where P 0 (X) is generally unknown but can be sampled from, e.g., via experimentally obtained statistics.
SI methods have been successfully applied to systems of biological relevance. Here we review two such cases involving 1) viral sequences [20][21][22][23] and 2) gene expression data [24]. A common challenge is the large number of features that are characteristic to the problem, namely, X = [x 1 , x 2 , ..., x q ] where q is large. The number of samples required to adequately sample a hyperspace grows exponentially with q, colloquially referred to as the "curse of dimensionality". [12] In such situations, where it is only possible to sample a small fraction of the possible feature combinations, the data will necessarily be sparse. The problem then is to formulate a model that can take such sparse sampling and make probability predictions for any X.
While statistically significant sampling of the full coordinate space, X, is challenging in high dimension, it is generally possible to obtain good statistics on quantities that depend on only small subset of coordinates-providing some constraints on P (X). For example, statistically significant empirical calculation of one-and two-point probability distributions, P is often possible. (For example, the radial distribution function, g(r), can be reliably computed from an MD simulation that samples only a minuscule fraction of possible particle arrangements.) Other quantities such as the mean x i P 0 (X) and co-variance x T i x j P 0 (X) may also be reliably estimated. However, enforcing P (X) to replicate lower-order correlations is an under-determined problem. A potentially infinite number of P (X) can reproduce such lower-order statistics. The challenge is to decide which of the innumerable possibilities for P (X) is "best".
Generally speaking, the best model choice given limited information is the one that is the least predictive (or least "informative") while reproducing the statistics that are known for P (X). [36] For example, two of the infinite probability distributions for a six-sided die that are known to yield an average face value of 3.5 are 1) a fair six-sided die (i.e., each face (f = 1 − 6) is equally likely, or P (f ) = 1/6 for all f ), and 2) a die that only turns up a 3 or a 4 in equal proportion (P (3) = P (4) = 1/2, P (f ) = 0 where f = 1, 2, 5, 6). Intuitively, the P (f ) for a fair die contains the least information possible-the distribution is maximally "smeared out" across possible states; by contrast, the latter distribution predicts that the face value of the roll will either be a 3 or a 4. For a more complex constraint though, the "best" P (X) may not so intuitive. Therefore, we require a measure quantifying the information content of a distribution P (X).
The maximum entropy principle states that the best model of the many that meet the desired constraints is the one that maximizes the Shannon entropy [36][37][38][39][40] where X is shorthand for a summation over all possible X. In practice, maximum entropy is typically carried out as a two-step protocol. First, Eqn. 3.1 is analytically maximized subject to the constraints via the technique of Lagrange multipliers. This provides a functional form for the distribution that is parametrized by the undetermined Lagrange multipliers, which are then solved for given the known constraints.

Viral Fitness Landscapes
Preventative vaccination represents the most realistic way of combating viruses such as Hepatitis C virus (HCV) and human immunodeficiency virus (HIV-I). However, many viruses have high sequence mutability, providing pathways to escape from vaccine-induced immunological pressure. Therefore, effective vaccines might need to account for the entirety of the vast combinatorial "sequence space" available to the virus. Unfortunately, relatively limited amounts of statistical data (from experimental observations of virus sequences) are available in comparison to the high dimensionality of the problem. Therefore, SI strategies based on maximizing the Shannon entropy are appropriate for this task. The "fitness" associated with a specific sequence is determined by its ability to replicate. While this quantity can be measured by in vitro studies, existing experimental data of this type are sparse. Instead, clinical databases that quantify the frequency of given sequences are used as a proxy for the fitness; these data span a much larger range of sequences. The underlying, physically reasonable, assumption is therefore that "fitter" sequences are more prevalent in the real world [21]. The goal of a vaccine is to induce immunological pressure on fit sequences, thus driving the viral population into unfit sequence space where it is minimally damaging to the infected host.
Both HCV and HIV-1 models followed from maximization of the Shannon entropy (Eqn. 3.1) subject to reproducing the empirical one-and two-point sequence distributions P where K is the combined set of all one-and two-point dimensionless "energy" functions, K (i) (s i ) and K (i,j) (s i , s j ), that need to be determined such that P (i) (s i ) = P (i) 0 (s i ) and P (i,j) (s i , s j ) = P (i,j) 0 (s i , s j ). Formally, this model is identical to the infinite range Potts spin model in statistical physics. [41] It would be possible to match higher-order correlations as well via a multi-body energy expansion; however, the sample size requirements for adequate statistical estimation of correlations increase exponentially with the dimensionality of the correlation quantity. Attempting to exactly reproduce higher-order correlations that have not been adequately sampled will generally result in over-fitting and therefore reduced transferability to data that is not included in the modeling step. Nonetheless, the above models based on one-and two-body correlations constructed for viruses were found to also reproduce three-body correlations, even though these were not explicitly considered for modeling. [22] The models were validated by a multifaceted analysis which tested the ability of the model to 1) infer sequence fitness, 2) identify viable escape mutations for the virus, and 3) agree with known effective immune responses that target specific regions of the viral sequence. With respect to the first component, one expects measurements of the in vitro fitness, F , to be related to the model via lnF (S|K) ∝ −E(S|K) if prevalence and fitness are proportional. As shown in Fig. 3a, this relation is found to be reasonably well approximated in HCV-an analogous study showed this to hold for HIV as well. Other more specific tests (such as identifying escape mutations) demonstrated similar success-validating the model and its underlying assumption of fitness and prevalence. In addition to providing an avenue to deal with the high dimensionality of the problem, the use of a spin model for viral sequences provides the added benefit of interpretability by way of analogy to actual statistical mechanical systems. In modeling HIV-I, a phase transition was observed that might be a consequence of viral "error catastrophe", a wellknown concept where accelerating the mutation rate of the virus to a certain tipping point eventually leads to viral extinction. [23] Above the transition, fatal error accumulation leads to an ensemble populated by many high energy (and therefore low fitness) strains. RNA viruses in particular are thought to have evolved to operate in close proximity to the error catastrophe, so as to be maximally mutable-enhancing the ability of viruses to escape the host's immune response-while remaining effective. [42] Purely theoretical models have been developed that display the error catastrophe; however, such models necessarily invoke uncontrolled assumptions about viral replication. [43] Empirical modeling of viral fitness as described above, on the other hand, is based upon clinical data that implicitly accounts for factors related to viral fitness and replication.
By defining a fictive temperature (T ), which corresponds to a re-scaling of the energy via E(S) → E(S)/T in Eqn. 3.2, a corresponding "heat capacity" was computed for the HIV-I fitness model. The sharp spike shown in Fig. 3b is suggestive of an underlying phase transition. This phase transition was interpreted as a signature of the error catastrophe because it signifies a change in the population preference for low versus high fitness sequences, as can be observed in Fig. 3c. Specifically, below the transition temperature (T ec = 1.20), the model virus exists as an ensemble of relatively few high fitness (i.e., low energy) strains, with the distribution narrowing and shifting to lowering energies as T decreases. Conversely, the model predicts a generally broader and much higher energy distribution above T ec . However, at T ec = 1.20, there is a co-existence between high and low fitness strains, as evidenced by the bimodality in the black curve in Fig. 3c. Of course, this co-existence between two "states" is also a well-known characteristic of a first-order phase change in systems that obey equilibrium statistical mechanics.
This research provides the first empirical HIV model to corroborate the limited experimental reports supporting the existence of an accessible error catastrophe. It has been proposed that altering the model via tuning T corresponds to varying the mutation rate: 1) T → ∞ corresponds to random copying, and 2) T → 0 corresponds to perfect, high-fidelity replication-a variable which could be exploited in the search for a cure. Interestingly, such effective temperatures may be reachable by mutation rate enhancing drug therapies, and some clinical trials for an HIV mutagen have already begun. [44][45][46]

Maximum entropy inference for complex systems
In the studies above, it might seem surprising that inclusion of only one-and two-body correlations is sufficient to accurately model the many-body relationship between sequence and viral fitness. However, recent work indicates that the complexity inherent to biological systems might actually contribute to the successes of pairwise models in approximating their behavior. [47] A prior study of simplified Ising p-spin models [48] indicates that the sufficiency of a pairwise model is tied to the entropy of the distribution being approximated. p-spin models describe a set of N q spins Q ≡ [q 1 , q 2 , ..., q Nq ] with possible values q i = {+1, −1} as follows where J p is shorthand for the set of all p-body spin couplings, {J ν 1 ,ν 2 ,...,νp }.
An ensemble of p = 3 and 4 models were randomly generated-possessing a range of entropies-by varying the number of spins (N q ), the number of non-zero p-body couplings (M J ), and the values of the couplings (J p ). These higher-order models were mapped onto pairwise spin models (p = 2), subject to the constraint that the two-body spin probabilities of the pairwise model matched the reduced two-body distributions obtained from the higherorder probability distribution. A schematic representation of this mapping for p = 3 is shown in Fig. 4a: 4). b) Randomly generated 3-body interactions networks composed of 11 ≤ N q ≤ 20 spins are mapped onto one-and two-body reduced models. c) Randomly generated 4-body interactions networks composed of 11 ≤ N q ≤ 22 spins are mapped onto a pairwise model. The data are segregated according to a measure of the network "interaction strength" α, which is proportional to the product of the standard deviation of the randomly generated spin couplings, s, and the (normalized) number of non-zero p-body couplings i.e., α = sM J /N q . Figure adapted from [47].
In order to evaluate how well the reduced distributions replicated the targets, a measure of similarity between two distributions is required. Various metrics exist for this purpose, many of which have seen renewed interest for their utility in rigorous force-field optimization and systematic coarse-graining. [49][50][51] For the p-spin study, the authors adopted the popular Kullback-Leibler divergence which is a measure of the degree to which P (X) overlaps with P 0 (X). [50] The KLdivergence is also known as the relative entropy in the field of bio-molecular coarse-graining. [50] Adopting the higher-order model as P 0 (X), Fig. 4b demonstrates that there are two regimes where the higher-order distributions are well-represented by pairwise models (i.e., D KL ≈ 0): when the entropy of the target is either very low or very high. For high entropy systems, this agreement is trivial: if there is relatively little order, there are not particularly stringent requirements on the (necessarily weak) couplings that can generate such a system. Even one-body, or independent, models work well when the entropy is sufficiently high, as can be seen in the right hand side of Fig. 4b. However, more surprising is the observation that very ordered systems (i.e., H → 0) are also well-represented by pairwise models. Furthermore, pairwise sufficiency was also found to be favored when the higher-order networks are densely interacting-as quantified by the parameter α. In Fig. 4c, the ensemble of spin models were first sorted by α of the target network and then plotted separately for different ranges in α. Fig. 4c shows a notable decrease in D KL as α increases.
Taken together, pairwise models are particularly well suited to represent higher-order networks when the latter has low entropy and is densely interacting. Low entropy systems necessarily possess correlations among constituent particles, spins, etc.; the lowest entropy spin system occurs when all spins move coherently (all up or all down), for instance. While such an arrangement may be obtained via many higher-order models, sufficiently strong positive pair interactions between all spins will achieve the same structure. Similarly, finite-sized correlated spin groups (i.e., clusters) can be obtained by employing strong pair interactions between the spins comprising each spin cluster while setting all other pair interactions to significantly smaller values. Interestingly, prior work on random constraint satisfaction anticipates the presence of such large correlated groups-amenable to representation via pair interactions-when the underlying distribution is densely interacting. [52,53] It has been suggested that biological systems can be regarded as evolutionarily derived solutions to complex constraint satisfaction problems [47,54]-ideal candidates for modeling with pairwise statistical models, according to the above analysis. Therefore, pairwise sufficiency may actually be characteristic to many biological systems-not just for the viral sequences reviewed earlier. Indeed, other successes in the modeling of biological processes with pairwise maximum entropy have been demonstrated. [40,[55][56][57][58] For example, such methods have been used to construct a model for gene expression levels under fluctuating cellular metabolic conditions. [24] DNA micro-arrays can measure thousands of gene expression levels at a time, resulting in a high-dimensional data set. Similar to the viral sequences above, pairwise entropy maximization provides an avenue to both circumvent the curse of dimensionality and extract pair interactions between genes.
This approach was applied to gene expression data of the well-studied Saccharomyces cerevisiae (i.e., Baker's yeast) under conditions which facilitate energy metabolic oscillations. [59] The resultant fluctuations in transcription levels yield an empirical ensemble of gene expression data. Micro-array data was encoded as an array of variables, G ≡ {g i }, where g i , is one of m continuous, scalar, gene expression levels. Constraints imposed on the maximum entropy model include the empirically observed one and two-body gene expression correlations, g i P 0 (G) and g i g j P 0 (G) , respectively; this yields a model analogous to Eqn. 3.3, where both p = 1 and 2 contributions are present and the spins are continuous.
The benefit of using SI for gene expression data as opposed to a more common correlationbased analysis (e.g., clustering) [60] is the construction of a network that quantifies the direct interactions between genes, possibly uncovering the driving force behind the observed correlations. Using correlation data alone only indicates which gene levels change together, but not what elements are driving those fluctuations. With the maximum entropy model though, it is easy to identify important governing genes, i.e., the nodes of high connectivity. For example, Fig. 5 has seven nodes with six or more connections, each of which corresponds to an important cellular protein for nutrient signaling. Furthermore, the model can uncover relationships between genes that are not part of the same cellular process. This can be seen visually in Fig. 5 by the many interactions between genes of different colors, where color indicates the cell process associated with a given gene. Therefore, while the protein expression patterns of S. cerevisiae are relatively well understood and characterized, this need not be the case for the maximum entropy model to discover gene interactions and "hubs" that govern cellular activity and system dynamics within the cell.

Statistical mechanical inverse design of structure and properties
The discovery of materials possessing a desired attribute is a long-standing design challenge that has been facilitated by computer simulation studies. Forward design-searching through parameter space to find favorable combinations that yield materials with some chosen characteristic-has benefited from recent advances in computational power. [61][62][63][64][65] However, combinatorial space increases exponentially with dimensionality; therefore, such approaches can still be hampered by the cost associated with exploration of vast swathes of parameter space as well as potential couplings between tunable parameters. Inverse design (ID) methodologies-those characterized by an optimization framework typically used to iteratively update the parameters in a systematic fashion-represent an alternative approach. Inverse methods of statistical mechanics have facilitated the discovery of conditions (particle interactions, external fields, etc.) that lead to assembly of novel architectures. [66][67][68][69][70] One challenge to ID is to delineate an effective scheme to update the parameters; SI techniques can be leveraged to construct an optimization framework for such purposes. For instance, the strategy outlined in Sect. 3 for modelling an unknown probability distribution function can be adapted for application to complex materials design problems by "mapping" a contrived probability distribution (P tgt (X))-which has been intentionally imbued with some desired property-onto a statistical mechanical model. Standard SI approaches can then be used to discover the optimal parameters needed to reproduce the property present in P tgt (X).
SI, therefore, can provide a prescription for performing updates to a probability distribution P (X) that is constrained by a statistical mechanical model. In the context of equilibrium molecular simulation, the clear choice for the probability of observing a given configuration X is the Boltzmann distribution The probability for the configuration X, where X can represent continuous particle coordinates, discrete spins, etc., can be tuned via the set of variables, θ ≡ {θ i } which parametrize the configuration energy U (X|θ). The focus of the remainder of this section is primarily on two recent advances in novel materials discovery, each pertaining to complementary aspects of materials design: structure and properties. Both strategies optimize a probability distribution toward the design goal, where the optimization is guided by SI.

Self-Assembly
Top down fabrication of materials with desired nanoscale structural characteristics remains a formidable technological challenge. [71] Spontaneous self-assembly of constituent nanoscale particles provides an alternative bottom up approach for realizing materials with such designer morphologies; however, careful tuning of inter-particle interactions is required. [61-65, 72, 73] A recent series of studies outlined a statistical mechanical ID strategy for the discovery of such interactions.
For a particle-based system, the relevant coordinate (X) is the set of D-dimensional Cartesian particle positions, R ≡ [r 1 , r 2 , ..., r N ]. Prior to performing SI, a contrived target probability distribution over configuration space, P tgt (R), that yields an ensemble of "desirable" configurations must be specified. [25][26][27][28][29] Simulations constrained by either manybody interactions or by use of multiple components that enforce the desired structure have been successfully employed for this purpose. Then the KL-divergence (Eqn. 3.4) between P (R|θ) and P tgt (R) is minimized (where P tgt (R) plays the role of P 0 (X) from Sect. 3). In Eqn. 3.4, the only term that depends on the tunable parameters θ is the second "overlap" contribution, lnP (X|θ) Ptgt(X) .
Maximizing lnP (R|θ) Ptgt(R) with respect to θ has an intuitive probabilistic interpretation under fairly nonrestrictive assumptions. In particular, it is equivalent to maximizing the probability for the parametric model, P (R|θ), to sample an independently identically distributed set of configurations [R 1 , R 2 , ..., R M ] drawn from P tgt (R) in the large sample limit (M → ∞)-formally known as the maximum-likelihood approach of SI. [11,12] Direct maximization of lnP (R|θ) Ptgt(R) is intractable due to the presence of Z(θ) in Eqn. 4.1; however, iterative updates to θ are easy to compute via natural gradient ascent. [25,26] Under the constraint of isotropic pair interactions, i.e., U (R|θ) ≡ N i<j=1 u(r i,j |θ) where r i,j ≡ |r i − r j |, gradient ascent optimization in the canonical ensemble corresponds to where g(r) is the radial distribution function, α is the "learning" rate which is empirically set to guarantee stability of the optimization, and the definition of Z(θ) from Eqn. 4.1 has been used. We note the analogy between the above approach to that outlined in Sect. 3; in both, pair correlations (here, g(r)) are used as a proxy for the full probability distribution. In a sense, MD configurations represent a familiar example of sparsely sampled data: there are infinite possibilities for particle arrangements, and only some finite subset are sampled in an MD simulation. The above approach has been successfully employed to engineer pair interactions that promote formation of clustered fluids, porous mesophase assemblies, and crystalline lattices; see Fig. 6. The strategy is sufficiently general to handle all of these disparate architectures, with the primary difference being the construction of appropriate target ensembles. In all cases, standard simulation tools were employed with constraints chosen to enforce desired structural motifs. For cluster fluids, a many-body constraint upon the radius of gyration of each cluster was employed in conjunction with inter-cluster repulsive interactions that effectively enforced a minimum distance between particles in different clusters. [25,29] In the case of a porous mesophase, a binary mixture of size-asymmetric hard-core-like particles was constructed, where the large particles created spherical zones of exclusion with respect to the smaller particles. [25,27] Finally, mutually non-interacting particles harmonically tethered to the desired lattice sites were used to emulate a crystalline lattice. [25,26] Though some of the above architectures had previously been observed in MD simulations with pairwise interactions, particularly for clustered phases, [74][75][76][77] no systematic framework existed for controlling the size of mesoscale structural features (i.e., clusters or the pores). By encoding the desired cluster or pore size in the target ensembles, interparticle interactions were discovered that reflected the prescribed mesophase dimension, allowing for the diameter of pores or the number of particles in a cluster to be tuned by modulating P tgt (R).
As can be visually gleaned from Fig. 6, the desired motifs (clusters, pores, and crystals) all require highly structured particle arrangements, intuitively indicating a target ensemble with relatively low entropy. As discussed in Sect. 3.2, such low entropy ensembles are generally amenable to pairwise models, and the successes outlined above seem to support this finding. [47] When there is flexibility in the creation of the target ensemble, a reasonable design principle might be to choose the lowest entropy ensemble that is practical. The systematic implementation of such a goal is an open question for future work.

Material Properties
While Sect. 4.1 outlines a robust methodology for designing structure, the relationship between structure and material properties is a complicated many-body physics problem. Therefore, a systematic, statistical mechanical framework for the direct optimization of targeted material properties is complementary to the structure-based ID approach outlined above.
One such ID scheme leverages the probabilistic micro-state configuration information afforded by a statistical mechanics driven approach. The key requirement of the method is a link between configuration, R, and a configuration-dependent FOM, F (R), which quantifies the user's goal. [30] For targeting a material property, F (R) could be the material property itself or a related metric. Derivation of a heuristic optimization strategy starts with an idealized case where P (R|θ) is infinitely flexible and can be additively updated by some δP id (R|θ) at each individual R. While the choice of δP id (R|θ) is nonunique, it is constrained by the requirement that normality of the probability is preserved, R δP id (R|θ) = 0, and that it generally increases (decreases) the probability for R that has improved (decreased) fitness. One simple choice that guarantees both is Constraining the probability distribution to conform to a realistic statistical mechanical model necessarily limits the possibilities for P (R|θ); i.e., the probability of observing a single configuration cannot generally be tuned independently of the probabilities associated with other configurations. Thus, exactly realizing Eqn. 4.3 is not generally feasible. Instead, one seeks the best approximation to the ideal update afforded by a change in parameters (δθ) via by minimizing the average (over P (R|θ)) sum of square differences between the ideal and realizable updates. Similar to Sect. 4.1, the updates to the parameters can be computed iteratively. The above methodology has been shown to be successful in a variety of contexts, including for a square-lattice Ising model where interaction parameters were tuned to maximize magnetization, optimization of a trap for a thermalized particle constrained to a 2D substrate, self-assembly of block copolymers on a chemically patterned substrate, and discovery of pairwise interactions that result in a six-bead chain spontaneously folding into an octahedron. For the last case, each bead, a coarse-grained representation of a polymer, was modeled as a hard sphere connected by fixed bond lengths between nearest neighbors. In-stead of optimizing for structure directly as in Sect. 4.1, the radius of gyration [R g (R)] of the chain was used to quantify the polymer shape. The qualitative progress of the optimization, progressing from an extended chain to a compact object, can be seen in Fig. 7a, where the index indicates the optimization step. More quantitatively, the metric R g /R oct − 1 (percent deviation from target, essentially an inverse fitness) decreased monotonically over the optimization; see Fig. 7b, where R g ≡ R g (R) P (R|θ) is the mean radius of gyration. After roughly 200 iterations, the mean deviation from an octahedral R g was only around 1%; it is visually apparent in Fig. 7a that the final steps in the optimization resulted in an octahedron.
The evolution of the couplings (λ) over the course of the optimization are shown in Fig. 7c. The inset schematic shown in Fig. 7c summarizes the optimized results as the interaction network where thicker lines indicate stronger couplings. Though symmetry is not explicitly enforced, the mirror symmetry present in both the bead model and the octahedron is manifest (to a good approximation) in the interactions, i.e., λ 1:5 ≈ λ 2:6 , etc. Furthermore, the resulting interactions are relatively intuitive; for instance, pairs of beads that are positioned co-linearly with the center of the octahedron, and therefore are farther apart, have the weakest couplings (λ 1:3 , λ 2:5 , λ 4:6 ).
For the above proof-of-concept example, the relationship between structure and property (R g ) is rather trivial, but this need not be the case. In principle, the method can be straightforwardly applied for any material property which uniquely depends on particle coordinates. For example, other structural metrics like the fractal dimension and average spacing between particles are feasible. More generally, quantities that use particle coordinates in intermediate theoretical calculations, such as electron hopping conductivity, [78] can be optimized.

Conclusions and outlook
Techniques of ML and SI are useful in a wide variety of physical chemistry applications in both experimental and computational contexts-a sampling of which has been provided in this review. For example, Bayesian optimization is straightforward to incorporate into experimental design loops to systematize the selection of new samples to prepare and measure. On the computational side, ML tools have been incorporated into MD simulations to accelerate QM force calculations, though other expensive computations for simulation could potentially be facilitated by ML as well. Also reviewed were applications of SI for modeling complex systems, with case studies taken from 1) biology in the form of viral sequences and gene expression data and 2) materials design of both structure and properties.
It is relatively early on in the adoption of ML and SI tools in the physical sciences. Undoubtedly many new applications remain to be discovered. However, one broad implication of further integrating ML and SI into research is the promise of increased opportunities for cooperation between experimental and computational research. For example, the viral sequence study used real clinical data to create an empirical computational model that may be translated into rational, real world, vaccine design. Similarly, Bayesian optimization can leverage theoretical frameworks to construct features to be input into a numerical optimiza-tion procedure that guides experimental design. In the future, perhaps the primary utility of ML and SI methodologies will be the erosion of the distinctions between experiment, computation, and theory, ultimately facilitating both collaboration and discovery.