Inferring Pairwise Interactions from Biological Data Using Maximum-Entropy Probability Models

Maximum entropy-based inference methods have been successfully used to infer direct interactions from biological datasets such as gene expression data or sequence ensembles. Here, we review undirected pairwise maximum-entropy probability models in two categories of data types, those with continuous and categorical random variables. As a concrete example, we present recently developed inference methods from the field of protein contact prediction and show that a basic set of assumptions leads to similar solution strategies for inferring the model parameters in both variable types. These parameters reflect interactive couplings between observables, which can be used to predict global properties of the biological system. Such methods are applicable to the important problems of protein 3-D structure prediction and association of gene–gene networks, and they enable potential applications to the analysis of gene alteration patterns and to protein design.


Introduction
Modern high-throughput techniques allow for the quantitative analysis of various components of the cell. This ability opens the door to analyzing and understanding complex interaction patterns of cellular regulation, organization, and evolution. In the last few years, undirected pairwise maximum-entropy probability models have been introduced to analyze biological data and have performed well, disentangling direct interactions from artifacts introduced by intermediates or spurious coupling effects. Their performance has been studied for diverse problems, such as gene network inference [1,2], analysis of neural populations [3,4], protein contact prediction [5][6][7][8], analysis of a text corpus [9], modeling of animal flocks [10], and prediction of multidrug effects [11]. Statistical inference methods using partial correlations in the context of graphical Gaussian models (GGMs) have led to similar results and provide a more intuitive understanding of direct versus indirect interactions by employing the concept of conditional independence [12,13].
Our goal here is to derive a unified framework for pairwise maximum-entropy probability models for continuous and categorical variables and to discuss some of the recent inference approaches presented in the field of protein contact prediction. The structure of the manuscript is as follows: (1) introduction and statement of the problem, (2) deriving the probabilistic model, (3) inference of interactions, (4) scoring functions for the pairwise interaction strengths, and (5) discussion of results, improvements and applications.
Better knowledge of these methods, along with links to existing implementations in terms of software packages, may be helpful to improve the quality of biological data analysis compared to standard correlation-based methods and increase our ability to make predictions of interactions that define the properties of a biological system. In the following, we highlight the power of inference methods based on the maximum-entropy assumption using two examples of biological problems: inferring networks from gene expression data and residue contacts in proteins from multiple sequence alignments. We compare solutions obtained using (1) correlation-based inference and (2) inference based on pairwise maximum-entropy probability models (or their incarnation in the continuous case, the multivariate Gaussian distribution).

Gene association networks
Pairwise associations between genes and proteins can be determined by a variety of data types, such as gene expression or protein abundance. Association between entities in these data types are commonly estimated by the sample Pearson correlation coefficient computed for each pair of variables x i and x j from the set of random variables x 1 ,. . ., x L . In particular, for M given samples in L measured variables, j À x j Þ denotes the (i, j)-element of the empirical covariance matrixĈ ¼ ðĈ ij Þ i;j¼1;...;L . The sample mean operator Á provides the empirical mean from the measured data and is defined as A simple way to characterize dependencies in data is to classify two variables as being dependent if the absolute value of their correlation coefficient is above a certain threshold (and independent otherwise) and then use those pairs to draw a so-called relevance network [14]. However, the Pearson correlation is a misleading measure for direct dependence as it only reflects the association between two variables while ignoring the influence of the remaining ones. Therefore, the relevance network approach is not suitable to deduce direct interactions from a dataset [15][16][17][18]. The partial correlation between two variables removes the variational effect due to the influence of the remaining variables (Cramér [19], p. 306). To illustrate this, let's take a simplified example with three random variables x A , x B , x C . Without loss of generality, we can scale each of these variables to zero-mean and unit-standard deviation by , which simplifies the correlation coefficient to r ij x i x j . The sample partial correlation coefficient of a three-variable system between x A and x B given x C is then defined as [19,20] The latter equivalence by Cramer's rule holds if the empirical covariance matrix, C ¼ ðĈ ij Þ i;j2fA;B;Cg , is invertible. Krumsiek et al. [21] studied the Pearson correlations and partial correlations in data generated by an in silico reaction system consisting of three components A, B, C with reactions between A and B, and B and C (Fig 1A). A graphical comparison of Pearson's correlations, r AB , r AC r BC , versus the corresponding partial correlations, r ABÁC , r ACÁB , r BCÁA , shows that variables A and C appear to be correlated when using Pearson's correlation as a dependency measure since both are highly correlated with variable B, which results in a false inferred reaction r AC . The strength of the incorrectly inferred interaction can be numerically large and therefore particularly misleading if there are multiple intermediate variables B [22]. The partial correlation analysis removes the effect of the mediating variable(s) B and correctly recovers the underlying interaction structure. This is always true for variables following a multivariate Gaussian distribution, but also seems to work empirically on realistic systems as Krumsiek et al. [21] have shown for more complex reaction structures than the example presented here.

Protein contact prediction
The idea that protein contacts can be extracted from the evolutionary family record was formulated and tested some time ago [23][24][25][26]. The principle used here is that slightly deleterious mutations are compensated during evolution by mutations of residues in contact in order to maintain the function and, by implication, the shape of the protein. Protein residues that are close in space in the folded protein are often mutated in a correlated manner. The main problem here is that one has to disentangle the directly co-evolving residues and remove transitive correlations from the large number of other co-variations in protein sequences that arise due to statistical noise or phylogenetic sampling bias in the sequence family. Interactions not internal to the protein are, for example, evolutionary constraints on residues involved in Reaction system reconstruction and protein contact prediction. Association results of correlationbased and maximum-entropy methods on biological data from an in silico reaction system (A) and protein contacts (B). (A) Analysis by Pearson's correlation yields interactions associating all three compounds A, B, and C, in contrast to the partial correlation approach which omits the "false" link between A and C. (Fig 1A  based on [21].) (B) Protein contact prediction for the human RAS protein using the correlation-based mutual information, MI, and the maximum-entropy based direct information, DI, (blue and red, respectively). The 150 highest scoring contacts from both methods are plotted on the protein contacts from experimentally determined structure in gray. (Fig 1B based  oligomerization, protein-protein, protein-substrate interactions [6,27,28]. In particular, the empirical single-site and pair frequency counts in residue i and in residues i and j for elements s, o of the 20-element amino acid alphabet plus gap, f i (s) and f ij (s, o), are extracted from a representative multiple sequence alignment under applied reweighting to account for biases due to undersampling. Correlated evolution in these positions was analyzed, e.g., by [29], by using the mutual information between residue i and j, Although results did show promise, an important improvement was made years later by using a maximum-entropy approach on the same setup [5][6][7]30]. In this framework, the direct information of residue i and j was introduced by replacing f ij in the mutual information by P dir ij , where P dir ij ðs; oÞ ¼ 1 Z ij expðe ij ðs; oÞ þh i ðsÞ þh j ðoÞÞ andh i ðsÞ;h j ðoÞ and Z ij are chosen such that P dir ij , which is based on a pairwise probability model of an amino acid sequence compatible with the iso-structural sequence family, is consistent with the single-site frequency counts. In an approximative solution, [6,7] determined the contact strength between the amino acids s and o in position i and j, respectively, by e ij ðs; oÞ ' ÀðC À1 ðs; oÞÞ ij : Here, (C −1 (s,o)) ij denotes the inverse element corresponding to C ij (s,o) f ij (s,o) − f i (s) f j (o) for amino acids s, o from a subset of 20 out of the 21 different states (the so-called gauge fixing, see below). The comparison of contact prediction results based on MI-and DI-score for the RAS human protein on top of the actual crystal structure shows a much more accurate prediction result when using the direct information instead of the mutual information ( Fig 1B).
The next section lays the foundation to deriving maximum-entropy models for the two data types: continuous, as used in the first example, and categorical, as used in the second one. Subsequently, we will present inference techniques to solve for their interaction parameters.

Deriving the Probabilistic Model
Ideally, one would like to use a probabilistic model that is, on the one hand, able to capture all orders of directed interactions of all observables at play and, on the other hand, correctly reproduces the observed and to-be-predicted frequencies. However, this would require a prohibitively large number of observed data points. For this reason, we restrict ourselves to probabilistic models with terms up to second order, which we derive for continuous, real-valued variables, and extend this framework to models with categorical variables that are suitable, for example, to treat sequence information in the next section.

Model formulation for continuous random variables
We model the occurrence of sets of events in a particular biological system by a multivariate probability distribution P(x) of L random variables x = (x 1 ,. . .,x L ) T 2R L that is, on the one hand, consistent with the mean and covariance obtained from M observed data values x 1 ,. . .,x M and, on the other hand, maximizing the information entropy, S, to obtain the simplest possible probability model consistent with the data. At this point, each of the data's variables x i is continuously distributed on real values. In a biological example, these data originate from gene expression studies and each variable x i corresponds to the normalized mRNA level of a gene measured in M samples. As an example, a recent pan-cancer study of The Cancer Genome Atlas (TCGA) provided mRNA levels from M = 3,299 patient tumor samples from 12 cancer types [31]. The problem can be large, e.g., in the case of a gene-gene association study one has L % 20,000 human genes.
The first constraint on the unknown probability distribution, P: R L !R !0 is that its integral normalizes to 1, ð which is a natural requirement on any probability distribution. Additionally, the first moment of variable x i is supposed to match the value of the corresponding sample mean over M measurements in each i = 1,. . ., L, where we define the n-th moment of the random variable x i distributed by the multivariate probability distribution P as hx n i i :¼ ð x PðxÞx n i dx. Analogously, the second moment of the variables x i and x j and its corresponding empirical expectation is supposed to be equal, for i, j = 1,. . ., L. Taken together, Eqs 4 and 5 constrain the distribution's covariance matrix to be coherent to the empirical covariance matrix. Finally, the probability distribution should maximize the information entropy, with the natural logarithm ln. A well-known analytical strategy to find functional extrema under equality constraints is the method of Lagrange multipliers [32], which converts a constrained optimization problem into an unconstrained one by means of the Lagrangian L. In our case, the probability distribution maximizing the entropy (Eq 6) subject to Eqs 3-5 is found as the stationary point of the Lagrangian L ¼ LðPðxÞ; a; β; γÞ [33,34], The real-valued Lagrange multipliers α, β = (β i ) i = 1,. . ., L and γ = (γ ij ) i,j = 1,. . ., L correspond to the constraints Eqs 3, 4, and 5, respectively. The maximizing probability distribution is then found by setting the functional derivative of L with respect to the unknown density P(x) to zero [33,35], g ij x i x j ¼ 0: Its solution is the pairwise maximum-entropy probability distribution, which is contained in the family of exponential probability distributions and assigns a nonnegative probability to any system configuration x = (x 1 ,. . .,x L ) T 2R L . For the second identity, we introduced the partition function as normalization constant, It can be shown by means of the information inequality that Eq 8 is the unique maximum-entropy distribution satisfying the constraints Eqs 3-5 (Cover and Thomas [35], p. 410). Note that α is fully determined for given β = (β i ) and γ = (γ ij ) by the normalization constraint Eq 3 and is therefore not a free parameter. The right-hand representation of Eq 8 is also referred to as Boltzmann distribution. The matrix of Lagrange multipliers γ = (γ ij ) has to have full rank in order to ensure a unique parametrization of P(x), otherwise, one can eliminate dependent constraints [33,36]. In addition, for the integrals in Eqs 3-6 to converge with respect to L-dimensional Lebesgue measure, we require γ to be negative definite, i.e., all of its eigenvalues to be negative or

Concept of entropy maximization
Shannon states in his seminal work that information and (information) entropy are linked: the more information is encoded in the system, the lower its entropy [37]. Jaynes introduced the entropy maximization principle, which selects for the probability distribution that is (1) in agreement with the measured constraints and (2) contains the least information about the probability distribution [38][39][40]. In particular, any unnecessary information would lower the entropy and, thus, introduce biases and allow overfitting. As demonstrated in the section above, the assumption of entropy maximization under first and second moment constraints results in an exponential model or Markov random field (in log-linear form) and many of the properties shown here can be generalized to this model class [41]. On the other hand, there is some analogy of entropy as introduced by Shannon to the thermodynamic notion of entropy. Here, the Second law of Thermodynamics states that each isolated system monotonically evolves in time towards a state of maximum entropy, the equilibrium. A thorough discussion of this analogy and its limitation in non-equilibrium systems is beyond the scope of this review, but can be found in [42,43]. Here, we exclusively use the notion entropy maximization as the principle of minimal information content in the probability model consistent with the data.

Categorical random variables
In the following section, we derive the pairwise maximum-entropy probability distribution on categorical variables. For jointly distributed categorical variables x = (x 1 ,. . .,x L ) T 2Ω L , each variable x i is defined on the finite set Ω = {s 1 ,. . ., s q } consisting of q elements. In the concrete example of modeling protein co-evolution, this set contains the 20 amino acids represented by a 20-letter alphabet from A standing for Alanine to Y for Tyrosine plus one gap element, then Ω = {A, C, D, E, F, G, H, I, K, L, M, N, P, Q, R, S, T, V, W, Y,−} and q = 21. Our goal is to extract co-evolving residue pairs from the evolutionary record of a given protein family. As input data, we use a so-called multiple sequence alignment, {x 1 ,. . ., x M } &Ω L×M , a collection of closely homologous protein sequences that is formatted such that it allows comparison of the evolution across each residue [44]. These alignments may stem from different hidden Markov model-derived resources, such as PFAM [45], hhblits [46], and Jackhmmer [47].
To formalize the derivation of the pairwise maximum-entropy probability distribution on categorical variables, we use the approach of [8,30,48] and replace, as depicted in Fig 2, each variable x i defined on categorical variables by an indicator function of the amino acid s 2 Ω, This embedding specifies a unique representation of any L-vector of categorical random variables, x, as a binary Lq-vector, x(σ) with a single non-zero entry in each binary q-subvector Inserting this embedding into the first and second moment constraints, corresponding to Eqs 3 and 4 in the continuous variable case, we find their embedded analogues, the single and pairwise marginal probability in positions i and j for amino acids s,o,2Ω and with the distribution's first moment in each random variable, hy i i ¼ X y PðyÞy i and y = (y 1 ,. . ., y Lq ) T 2R Lq . The analogue of the covariance matrix then becomes a symmetric Lq × Lq matrix of connected correlations whose entries C ij (s,o) = P ij (s,o) − P i (s) P j (o) characterize the dependencies between pairs of variables. In the same way, the sample means translate to the single-site and pair frequency counts over m = 1,. . .,M data vectors The pairwise maximum-entropy probability distribution in categorical variables has to fulfill the normalization constraint, X x PðxÞ ¼ Furthermore, the single and pair constraints, the analogues of Eqs 3 and 4, enforce the resulting probability distribution to be compatible with the measured single and pair frequency counts, for each i, j = 1,. . ., L and amino acids s,o2Ω. As before, we require the probability distribution to maximize the information entropy, The corresponding Lagrangian, L ¼ LðPðxðσÞÞ; a; βðσÞ; γðσ; ωÞÞ, has the functional form, i;j¼1;...;L , respectively. The Lagrangian's stationary point, found as the solution of @L @PðxðσÞÞ ¼ 0, determines the pairwise maximum-entropy probability distribution in categorical variables [30,49], This is the 21-state maximum-entropy probability distribution as presented by [5][6][7].

Gauge fixing
In contrast to the continuous variable case in which the number of constraints naturally matches the number of unknown parameters, the case of categorical variables has dependencies due to 1 ¼ X s2O P i ðsÞ for each i = 1,. . ., L and P i ðsÞ ¼ X o2O P ij ðs; oÞ for each i, j = 1,. . ., L and s2Ω. This results in at most LðLÀ1Þ 2 ðq À 1Þ 2 þ Lðq À 1Þ independent constraints compared to LðLÀ1Þ 2 q 2 þ Lq free parameters to be estimated. To ensure the uniqueness of the inferred parameters in defining the Hamiltonian, Hðx 1 ; . . . ; x L Þ ¼ À X i<j e ij ðx i ; x j Þ À X i h i ðx i Þ, and, by implication, the probability distribution, one has to reduce the number of independent parameters such that these match the number of independent constraints. For this purpose, so-called gauge fixing [5] has been proposed, which can be realized in different ways. For example, the authors of [6,7] set the parameters corresponding to the last amino acid in the alphabet, s q , to zero, i.e., e ij (s q ,Á) = e ij (Á, s q ) = 0 and h i (s q ) = 0 for 1 i < j L, resulting in rows and columns of zeros at the end of each q Â q-block of the Lq × Lq coupling matrix. Alternatively, the authors of [5] introduce a zero-sum gauge,

Network interpretation
The derived pairwise maximum-entropy distributions in Eqs 13 or 12 and 8 specify an undirected graphical model or Markov random field [34,41]. In particular, a graphical model represents a probability distribution in terms of a graph that consists of a node and an edge set. Edges characterize the dependence structure between nodes and a missing edge then corresponds to conditional independence given the remaining random variables. For continuous, real-valued variables, the maximum-entropy distribution with first and second moment constraints is multivariate Gaussian, which will be demonstrated in the next section. Its dependency structure is represented by a graphical Gaussian model (GGM) in which a missing edge, γ ij = 0, corresponds to conditional independence between the random variables x i and x j (given the remaining ones), and is further specified by a zero entry in the corresponding inverse covariance matrix, (C −1 ) ij = 0.
In the next section, we describe how the dependency structure of the graph is inferred.

Inference of Interactions
Up to this point, the functional form of the maximum-entropy probability distribution is specified, but not its determining parameters. For categorical variables with dimension L > 1, there is typically no closed-form solution. In the following section, we present several inference methods to estimate these parameters that have recently been used in the context of protein contact prediction. Those are (1) for continuous variables, the exact closed-form solution which approximates the mean-field result for categorical variables, and (2) three inference methods for categorical variables based on the maximum-likelihood methodology: the stochastic maximum likelihood, the approximation by pseudo-likelihood maximization, and finally, the sparse maximum-likelihood solution.

Closed-Form Solution for Continuous Variables
The simplest approach to extract the unknown Lagrange multipliers α, β = (β i ), and γ = (γ ij ) from P(x) exactly is to use basic integration properties of the continuous random variables x i in the constraints Eqs 3-5. For this purpose, we rewrite the exponent of the pairwise maximumentropy probability distribution Eq 8, ðx Àγ À1 βÞ Tγ ðx Àγ À1 βÞ ; where we use the replacementγ :¼ À2γ and requireγ to be positive definite (which is equivalent to γ being negative definite), i.e., x Tγ x > 0 for any x 6 ¼ 0, which makes its inverseγ À1 ¼ À 1 2 γ À1 well-defined. As already discussed, this is a sufficient condition on the integrals in Eqs 3-6 to be finite. For notational convenience, we define the shifted variable z ¼ ðz 1 ; . . . ; z L Þ T :¼ x Àγ À1 β or x i ¼ z i þ X L j¼1 ðγ À1 Þ ij b j and accordingly, the maximum-entropy probability distribution becomes with the normalization constantZ ¼ exp 1 À a À 1 2 β Tγ À1 β À Á . The normalization condition Eq 3 in the new variable is, and the linear shift does not affect the integral when integrated over R L yielding for the normalization constant,Z ¼ ð z e À 1 2 z Tγ z dz. Furthermore, the first-order constraint Eq 4 becomes for each i = 1,. . ., L, and we used the point symmetry of the integrand then, Analogously, we find for the second moment, determining the correlations for each index pair i, j = 1,. . ., L, where we use again the point symmetry and the result on the normalization constraint. Based on this, the covariance is found as, Finally, the term hz i z j i is solved using a spectral decomposition of the symmetric and positive-definite matrixγ as sum over products of its eigenvectors v 1 ,. . .,v L and real-valued and pos- The eigenvectors form a basis of R L and assign new coordinates, y 1 ,. . .,y L , to z ¼ X L k¼1 y k v k , which allows writing of the exponent hz i z j i as z Tγ z ¼ X L k¼1 l k y 2 k . The covariance between x i and x j then reads as (Bishop [52], p. 83) Taken together, the Lagrange multipliers β and γ are specified in terms of the mean, hxi, and the inverse covariance matrix (also known as the precision or concentration matrix), C −1 , As a consequence, the real-valued maximum-entropy distribution Eq 14 for given first and second moments is found as the multivariate Gaussian distribution, which is determined by the mean hxi and the covariance matrix C, Pðx; hxi; CÞ ¼ ð2pÞ ÀL=2 detðCÞ À1=2 exp À 1 2 ðx À hxiÞ T C À1 ðx À hxiÞ ð17Þ and we refer to [52] for the derivation of the normalization factor. The initial requirement of γ ¼ À2γ to be positive definite results in a positive-definite covariance matrix C, a necessary condition for the Gaussian density to be well defined. In summary, the multivariate Gaussian distribution maximizes the entropy among all probability distributions of continuous variables with specified first and second moments. The pair interaction strength is now evaluated by the already introduced partial correlation coefficient between x i and x j given the remaining variables {x r } r2{1,. . ., L}\{i,j} , r ijÁf1;...;Lgnfi;jg g ij ffiffiffiffiffiffiffiffi

Data integration
In biological datasets as used to study gene association, the number of measurements, M, is typically smaller than the number of observables, L, i.e., M < L in our terminology. Consequently, the empirical covariance matrix,Ĉ ¼ 1 M X M m¼1 ðx m À xÞðx m À xÞ T , will in these cases always be rank-deficient (and, thus, not invertible) since its rank can exceed neither the number of variables, L, nor the number of measurements, M. Moreover, even in cases when M ! L, the empirical covariance matrix may become non-invertible or badly conditioned (i.e., close to singular) due to dependencies in the data. However, for variables following a multivariate Gaussian distribution, one can access the elements of its inverse by maximizing the penalized Gaussian loglikelihood, which results in the following estimate of the inverse covariance matrix, C À1 % C À1 d;l , C À1 d;l ¼ arg max with penalty parameter λ ! 0 and kYk d d ¼ If λ = 0, we obtain the maximum-likelihood estimate, for δ = 1 and λ > 0 the ℓ 1 -regularized (sparse) maximum-likelihood solution that selects for sparsity [53,54], and for δ = 2 and λ > 0 the ℓ 2 -regularized maximum-likelihood solution that favors small absolute values in the entries of the selected inverse covariance matrix [55]. For δ = 1 and λ > 0, the method is called LASSO, for δ = 2 and λ > 0, ridge regression. Alternatively, regularization can be directly applied to the covariance matrix, e.g., by shrinkage [17,56].

Solution for categorical variables
An ad hoc ansatz to extract the pairwise parameters in the categorical variables case (12) is to extend the binary variable xðσÞ ¼ ðx i ðs k ÞÞ iÁk 2 f0; 1g LðqÀ1Þ to a continuous one, y = (y j ) j 2R L(q−1) , and replace the sums in the distribution and the moments hÁi by integrals. The extended binary maximum-entropy distribution Eq 12 is then approximated by the Lq-dimensional multivariate Gaussian with inherited analogues of the mean hyi ¼ ðf i ðs k ÞÞ iÁk 2 R LðqÀ1Þ and the empirical covariance matrixĈðσ; ωÞ ¼ ðĈ ij ðs k ; s l ÞÞ i;j;k;l 2 R LðqÀ1ÞÂLðqÀ1Þ whose elementsĈ ij ðs; oÞ ¼ f ij ðs; oÞ À f i ðsÞf j ðoÞ are characterizing the pairwise dependency structure. The gauge fixing results in setting the preassigned entries referring to the last amino acid in the mean vector and the covariance matrix to zero, which reduces the model's dimension from Lq to L(q−1); otherwise the unregularized covariance matrix would always be non-invertible. Typically, the single and pair frequency counts are reweighted and regularized by pseudocounts (see section "Sequence data preprocessing") to additionally ensure thatĈðσ; ωÞ is invertible. Final application of the closed-form solution for continuous variables Eq 16 to the extended binary variables for C À1 ðσ; ωÞ %Ĉ À1 ðσ; ωÞ yields the so-called mean-field (MF) approximation [48], g MF ij ðs; oÞ ¼ À for amino acids s,o2Ω and with restriction to residues i < j in the latter identity. The same solution has been obtained by [6,7] using a perturbation ansatz to solve the q-state Potts model termed (mean-field) Direct Coupling Analysis (DCA or mfDCA). In Ising models, this result is also known as naïve mean-field approximation [57][58][59].
The following section is dedicated to maximum likelihood-based inference approaches, which have been presented in the field of protein contact prediction.

Maximum-Likelihood Inference
A well-known approach to estimate the parameters of a model is maximum-likelihood inference. The likelihood is a scalar measure of how likely the model parameters are, given the observed data (Mackay [34], p. 29), and the maximum-likelihood solution denotes the parameter set maximizing the likelihood function. For Markov random fields, the maximum-likelihood solution is consistent, i.e., recovers the true model parameters in the limit of infinite data (Koller and Friedman [32], p. 949). In particular, for a pairwise model with parameters hðσÞ ¼ The estimates of the model parameters are then obtained as the maximizer of l or, using the monotonicity of the logarithm, the minimizer of ln l, fh ML ðσÞ; e ML ðσ; ωÞg ¼ arg max hðsÞ;eðs;oÞ lðhðσÞ; eðσ; ωÞÞ arg min hðsÞ;eðs;oÞ À ln lðhðσÞ; eðσ; ωÞÞ: When we specify the maximum-entropy distribution Eq 13 as model distribution, the thenconcave loglikelihood [32] becomes ln lðhðσÞ; eðσ; ωÞÞ ¼ ln Pðx m ; hðσÞ; eðσ; ωÞÞ The maximum-likelihood solution is found by taking the derivatives of Eq 22 with respect to the model parameters h i (s) and e ij (s,o) and setting to zero, @ @h i ðsÞ ln l ¼ ÀM @ @h i ðsÞ ln Z j fhðσÞ;eðσ;ωÞg À f i ðsÞ ! ¼ 0; @ @e ij ðs; oÞ ln l ¼ ÀM @ @e ij ðs; oÞ ln Z j fhðσÞ;eðσ;ωÞg À f ij ðs; oÞ The partial derivatives of the partition function, in residues i = 1,. . ., L and i,j = 1,. . ., L, respectively, and for amino acids s,o2Ω. In other words, matching the moments of the pairwise maximum-entropy probability distribution to the given data is equivalent to maximum-likelihood fitting of an exponential family [34,60]. Although the maximum-likelihood solution is globally optimal for the pairwise maximumentropy probability model, based on the concavity of ln l, the resulting distribution is not necessarily unique, due to dependencies in the input data (Koller and Friedman [32], p. 948). To remove these equivalent optima and select for a unique representation, one needs to introduce further constraints by, for example, gauge fixing or regularization.
Each factor is of the following analytical form, which only depends on the unknown parameters (e ij (s,o)) i6 ¼r,j6 ¼r and (h i (s)) i6 ¼r and makes the computation of the pseudo-likelihood tractable. Note, we treat e ij (s,o) = e ji (o,s) and e ii (Á,Á) = 0. By this approximation, the loglikelihood Eq 21 becomes the pseudo-loglikelihood, ln l PL ðhðσÞ; eðσ; ωÞÞ :¼ ln Pðx r ¼ x m r jx nr ¼ x m nr ; hðσÞ; eðσ; ωÞÞ: In the final formulation of the pseudo-likelihood maximization (PLM) problem, an ℓ 2 -regularizer is added to select for small absolute values of the inferred parameters, fh PLM ðσÞ; e PLM ðσ; ωÞg ¼ arg min hðσÞ;eðσ;ωÞ fÀln l PL ðhðσÞ; eðσ; ωÞÞ þ l h khðσÞk 2 2 þl e keðσ; ωÞk 2 2 g; where λ h , λ e > 0 adjust the complexity of problem and are selected in a consistent manner across different protein families to avoid overfitting. This approach has been presented (with scaling of the pseudo-loglikelihood by 1 M eff w m to include sequence weighting, see section "Sequence data preprocessing") by [51] under the name plmDCA (PseudoLikelihood Maximization Direct Coupling Analysis) and has shown performance improvements compared to the mean-field approximation Eq 20. Another inference method based on the pseudolikelihood maximization but including prior knowledge in terms of secondary structure and information on pairs likely to be in contact is Gremlin (Generative REgularized ModeLs of proteINs) [65][66][67].

Sparse maximum likelihood
Similar to the derivation of the mean-field result (20), Jones et al. [8] approximated Eq 12 by a multivariate Gaussian and accessed the elements of the inverse covariance matrix by a maximum-likelihood inference under sparsity constraint [54,68,69]. The corresponding method has been called Psicov (Protein Sparse Inverse COVariance). The validity of this approach to solve the sparse maximum-likelihood problem in binary systems such as Ising models has been demonstrated by [69], followed by consistency studies [70]. In particular, the Psicov method infers the sparse maximum-likelihood estimate of the inverse covariance matrix Eq 19 for δ = 1 using the analogue of the empirical covariance matrix derived from the observed amino acid frequencies,Ĉðσ; ωÞ. Its elementsĈ ij ðs; oÞ ¼ f ij ðs; oÞ À f i ðsÞf j ðoÞ, the empirical connected correlations, are preprocessed by reweighting and regularized by pseudocounts and shrinkage. Regularized loglikelihood maximization Eq 19 selects a unique representation of the model, i.e., no additional gauge fixing is required. Using identity Eq 16 on the elements of the sparse maximum-likelihood (SML) estimate of the inverse covariance, C À1 1;l ðσ; ωÞ, yields the estimates for the Lagrange multipliers, g SML ij ðs; oÞ ¼ À

Sequence data preprocessing
The study of residue-residue co-evolution is based on data from multiple sequence alignments, which represent sampling from the evolutionary record of a protein family. Multiple sequence alignments from currently existing sequence databases do not evenly represent the space of evolved sequences as they are subject to acquisition bias towards available species of interest. To account for uneven representation, sequence reweighting has been introduced to lower the contributions of highly similar sequences and assign higher weight to unique ones (see Durbin et al. [44], p. 124 ff.). In particular, the weight of the m-th sequence, w m : = 1/k m , in the alignment {x 1 ,. . .,x M }, can be chosen to be the inverse of k m :¼ the number of sequences x m shares more than θ Á 100% of its residues with. Here, θ denotes a similarity threshold and is typically chosen as 0.7 θ 0.9, 1(a,b) = 1 if a = b and 1(a,b) = 0, otherwise, and H is the step function with H(y) = 0 if y < 0 and H(y) = 1, otherwise. This also provides us with an estimate of the effective number of sequences in the alignment, Additionally, pseudocount regularization withl > 0 is used to deal with finite sampling bias and to account for underrepresentation [5][6][7][8]44,48], resulting in zero entries inĈðσ; ωÞ, for instance, if a certain amino acid pair is never observed. The use of pseudocounts is equivalent to a maximum a posteriori (MAP) estimate under a specific inverse Wishart prior on the covariance matrix [48]. Both preprocessing steps combined yield the reweighted single and pair frequency counts, in residues i,j = 1,. . ., L and for amino acids s,o2Ω. Ideally for maximum-likelihood inference, the random variables are assumed to be independent and identically distributed. However, this is typically violated in realistic sequence data due to phylogenetic and sequencing bias, and the reweighting presented here does not necessarily solve this problem.

Scoring Functions for the Pairwise Interaction Strengths
For pairwise maximum-entropy models of continuous variables, the natural scoring function for the interaction strength between two variables x i and x j , given the inferred inverse covariance matrix, is the partial correlation Eq 18. However, for categorical variables, the situation is more complicated, and there are several alternative choices of scoring functions. Requirements on the scoring function are that it has to account for the chosen gauge and, in the case of protein contact prediction, evaluate the coupling strength between two residues i and j summarized across all possible q 2 amino acids pairs. The highest scoring residue pair is, for instance, used to predict the 3-D structure of the protein of interest. For this purpose, the direct information, defined as the mutual information applied to P dir ij ðs; oÞ ¼ : However, this expression is not gauge-invariant [5]. In this context, the notation with e ij (s, o), which refers to indices restricted to i < j, is extended and treated such that e ij (s,o) = e ji (o, s) and e ij (Á,Á) = 0; then ||e ij || F = ||e ji || F and ||e ii || F = 0. In order to correct for phylogenetic biases in the identification of co-evolved residues, Dunn et al. [27] introduced the average product correction (APC). It has been originally used in combination with the mutual information but was recently combined with the ℓ 1 -norm [8] and the Frobenius/ℓ 2 -norm [51] and is derived from the averages over rows and columns of the corresponding norm of the matrix of the e ij parameters. In this formulation, the pair scoring function is for e ij -parameters fixed by zero-sum gauge and with the means over the non-zero elements in row, column and full matrix, ke iÁ k F : ke ij k F and ke ÁÁ k F :¼ 1 LðLÀ1Þ X L i;j¼1 ke ij k F , respectively. Alternatively, the average product-corrected ℓ 1 -norm applied to the 20×20-submatrices of the estimated inverse covariance matrix, in which contributions from gaps are ignored, has been introduced by the authors of [8] as the Psicov-score. Using the average product correction, the authors of [51] showed for interaction parameters inferred by the mean-field approximation that scoring with the average product-corrected Frobenius norm increased the precision of the predicted contacts compared to scoring with the DI-score. The practical consequence of the choice of scoring method depends on the dataset and the parameter inference method.

Discussion of Results, Improvements, and Applications
Maximum entropy-based inference methods can help in estimating interactions underlying biological data. This class of models, combined with suitable methods for inferring their numerical parameters, has been shown to reveal-to a reasonable approximation-the direct interactions in many biological applications, such as gene expression or protein residue-residue coevolution studies. In this review, we have presented maximum-entropy models for the continuous and categorical random variable case. Both approaches can be integrated into a framework, which allows the use of solutions obtained for continuous variables as approximations for the categorical random variable case (Fig 3).
The validity and precision of the available maximum-entropy methods could be improved to yield more biologically insightful results in several ways. Advanced approximation methods derived from Ising model approaches [59,71] are possible extensions for efficient inference. Moreover, additional terms beyond pair interactions can be included in models of continuous and discrete random variables [1,33,59]. However, higher-order models demand more data, which is a major bottleneck for their application to biological problems. In the case of protein contact prediction, this could be resolved by getting more sequences, which are being obtained as the result of extraordinary advances in sequencing technology. The quality of existing methods can be improved by careful refinement of sequence alignments in terms of cutoffs and gaps or by attaching optimized weights to each of the data sequences. Alternatively, one could try to improve the existing model frameworks by accounting for phylogenetic progression [27,49,72] and finite sampling biases.
The advancement of inference methods for biological datasets could help solve many interesting biological problems, such as protein design or the analysis of multi-gene effects in relating variants to phenotypic changes as well as multi-genic traits [73,74]. The methods presented here could help reduce the parameter space of genome-wide association studies to first approximation. In particular, we envision the following applications: (1) in the disease context, coevolution studies of oncogenic events, for example copy number alterations, mutations, fusions Scheme of pairwise maximum-entropy probability models. The maximum-entropy probability distribution with pairwise constraints for continuous random variables is the multivariate Gaussian distribution (left column). For the maximum-entropy probability distribution in the categorical variable case (right column), various approximative solutions exist, e.g., the mean-field, the sparse maximum-likelihood, and the pseudolikelihood maximization solution. The mean-field and the sparse maximum-likelihood result can be derived from the Gaussian approximation of binarized categorical variables (thin arrow). Pair scoring functions for the continuous case are the partial correlations (left column). For the categorical variable case, the direct information, the Frobenius norm, and the average product-corrected Frobenius norm are used to score pair couplings from the inferred parameters (right column). and alternative splicing, can be used to derive direct co-evolution signatures of cancer from available data, such as The Cancer Genome Atlas (TCGA); (2) de novo design of protein sequences as, for example, described in [65,75] for the WW domain using design rules based on the evolutionary information extracted from the multiple sequence alignment; and (3) develop quantitative models of protein fitness computed from sequence information.
In general, in a complex biological system, it is often useful for descriptive and predictive purposes to derive the interactions that define the properties of the system. With the methods presented here and available software (Table 1), our goal is not only to describe how to infer these interactions but also to highlight tools for the prediction and redesign of properties of biological systems. pseudolikelihood maximization plmDCA APC-FN-score [78][79][80] pseudolikelihood maximization Gremlin Gremlin-score [81] sparse maximum-likelihood Psicov Psicov-score [82] continuous sparse maximum-likelihood glasso partial correlations [83] ℓ 2 -regularized maximum-likelihood scout partial correlations [84] shrinkage corpcor, GeneNet partial correlations [85,86] doi:10.1371/journal.pcbi.1004182.t001