Monte Carlo simulation of a statistical mechanical model of multiple protein sequence alignment

TL;DR

This paper introduces a Monte Carlo simulation method for a lattice gas model of multiple protein sequence alignment, enabling detailed exploration of sequence space with long-range interactions and insertions, and reveals a two-state transition in the SH3 domain family.

Contribution

It presents a novel grand canonical Monte Carlo algorithm for the lattice gas model of protein sequences, incorporating long-range interactions and variable insertions, and compares it with mean-field approximations.

Findings

Identification of a two-state transition in the SH3 domain family

Demonstration of the inaccuracy of mean-field approximation for the model

Preliminary results across various temperatures and chemical potentials

Abstract

A grand canonical Monte Carlo (MC) algorithm is presented for studying the lattice gas model (LGM) of multiple protein sequence alignment, which coherently combines long-range interactions and variable-length insertions. MC simulations are used for both parameter optimization of the model and production runs to explore the sequence subspace around a given protein family. In this Note, I describe the details of the MC algorithm as well as some preliminary results of MC simulations with various temperatures and chemical potentials, and compare them with the mean-field approximation. The existence of a two-state transition in the sequence space is suggested for the SH3 domain family, and inappropriateness of the mean-field approximation for the LGM is demonstrated.

Figures

Equations

$$E(\mathbf{X})=-\sum _{k=1}^{L_{\mathbf{X}}-1}J(X_k,X_{k+1})-\sum _{(s(X_k),s(X_l))\in \mathcal{T}}K(X_k,X_l)-\sum _{k=1}^{L_{\mathbf{X}}}\mu (X_k)$$

$$P(\mathbf{X})=\frac{\exp [-E(\mathbf{X})/T]}{\Xi [T]}$$

$$\Xi [T]=\sum _{\mathbf{X}}\exp [-E(\mathbf{X})/T].$$

$$\Omega [T]=-T\log \Xi [T].$$

$$n_Y(\mathbf{X})=\sum _{k=1}^{L_{\mathbf{X}}}{\delta }_{Y,X_k}$$

$$n_{Y,Y^{\prime }}^b(\mathbf{X})=\sum _{k=1}^{L_{\mathbf{X}}-1}{\delta }_{Y,X_k}{\delta }_{Y^{\prime },X_{k+1}}.$$

$$\sum _{a\in {\mathcal{A}}_{O_i}}{n}_{(O_i,a)}(\mathbf{X})$$

$$\sum _{b\in {\mathcal{A}}_{O_{i+1}}}{n}_{(S,a),(O_{i+1},b)}^b(\mathbf{X})+\sum _{b\in {\mathcal{A}}_{I_i}}{n}_{(S,a),(I_i,b)}^b(\mathbf{X})$$

$$\sum _{a\in {\mathcal{A}}_{O_i}}{n}_{(O_i,a),(S^{\prime },b)}^b(\mathbf{X})+\sum _{a\in {\mathcal{A}}_{I_i}}{n}_{(I_i,a),(S^{\prime },b)}^b(\mathbf{X})$$

$$\sum _{b\in {\mathcal{A}}_{O_j}}{n}_{(O_i,a)}(\mathbf{X})n_{(O_j,b)}(\mathbf{X})$$

$$\sum _{a\in {\mathcal{A}}_{O_i}}{n}_{(O_i,a)}(\mathbf{X})n_{(O_j,b)}(\mathbf{X})$$

$$E(\mathbf{X})$$

$$\min \{1,\exp [-(E_b-E_a)/T]\}.$$

$$\min \left\{1,\frac{p_{l+1}^{-}{q}_{21}}{p_l^{+}{q}_{20}}\exp [-(E_{l+1}-E_l)/T]\right\}$$

$$\min \left\{1,\frac{p_{l-1}^{+}{q}_{20}}{p_l^{-}{q}_{21}}\exp [-(E_{l-1}-E_l)/T]\right\}$$

$${⟨n_Y⟩}_{\mathrm{obs}}$$

$${⟨n_{Y{Y}^{\prime }}^b⟩}_{\mathrm{obs}}$$

$${⟨n_Y{n}_{Y^{\prime }}⟩}_{\mathrm{obs}}$$

$${⟨Q⟩}_{\mathrm{sim}}=\frac{1}{M}\sum _{t=1}^{M}Q({\mathbf{X}}^t).$$

$${⟨E⟩}_{\mathrm{obs}}=-\sum _{(Y,Y^{\prime })}^{b.p.}J(Y,Y^{\prime }){⟨n_{Y,Y^{\prime }}^b⟩}_{\mathrm{obs}}-\sum _{(Y,Y^{\prime })}^{n.b.p.}K(Y,Y^{\prime }){⟨n_Y{n}_{Y^{\prime }}⟩}_{\mathrm{obs}}.$$

$$F[T]=-\underset{J,K}{\max }\left({⟨E⟩}_{\mathrm{obs}}-\Omega [T]\right)$$

$${⟨Q⟩}_{\mathrm{obs}}$$

$$⟨Q⟩$$

$$R^{(\nu +1)}$$

$$P^{(\nu +1)}$$

$${\tilde{K}}_{(O_i,a)}=\sum _{(O_j,b)}K((O_i,a),(O_j,b)){⟨n_{(O_j,b)}⟩}_{\mathrm{mf}}$$

$${⟨n_{(O_j,b)}⟩}_{\mathrm{mf}}=\sum _{\mathbf{X}}{n}_{(O_j,b)}(\mathbf{X})P_{\mathrm{mf}}(\mathbf{X})$$

$$D_{O_i}=\sum _a⟨n_{(O_i,a)}⟩\ln \left[⟨n_{(O_i,a)}⟩/{⟨n_{(O_i,a)}⟩}_{\mathrm{obs}}\right]$$

$${\Delta }_{I_i}=\sqrt{\sum _a(⟨n_{(I_i,a)}⟩-{⟨n_{(I_i,a)}⟩}_{\mathrm{obs}}{)}^2}.$$

Full text

I Introduction

Thanks to the massive genome sequencing, we have a great number of known amino acid sequences at our hands. Exploiting the wealth of sequence data, recent advances in biological sequence analysis made it possible to reliably extract the “direct” couplings between residues that are separated along the sequence Morcos et al. (2011); Jones et al. (2012); Taylor et al. (2012); Miyazawa (2013); Levy et al. (2017), and thereby to accurately predict three-dimensional (3D) structures Marks et al. (2011); Ovchinnikov et al. (2017) as well as mutation effects Hopf et al. (2017).

Based on these developments, I have previously proposed a lattice gas model (LGM) of multiple protein sequence alignment (MSA) which incorporates direct couplings and variable-length insertions in a coherent manner Kinjo (2016). In that work, I tried to use a mean-field approximation for treating the long-range direct couplings so that the partition function can be computed efficiently. However, it was found that the quasi-1-dimensional model structure was somehow incompatible with long-range interactions so that the original mean-field approximation without the diagonal terms of the direct coupling matrix failed to converge to correct solutions. Thus, I resorted to the Gaussian approximation by including the diagonal terms of the “long-range” interaction matrix. The Gaussian approximation makes the system essentially harmonic, and hence, by construction, it does not exhibit some interesting phenomena such as phase transitions. In a preliminary study, I also tested the pseudolikelihood method Balakrishnan et al. (2011); Ekeberg et al. (2013) for obtaining the direct couplings, but again, failed to obtain stable solutions (data not shown). Therefore, for the LGM of protein families, it appears necessary to drop any approximations (at least those known to the author). In principle, this can be done by performing Monte Carlo (MC) simulations, which was indeed the original approach to the problem Lapedes et al. (1999), and more recently employed by Sutto et al.Sutto et al. (2015). In these studies, only alignments with fixed lengths were treated so that the standard canonical Metropolis sampling was sufficient. In the LGM, however, the alignment length is variable, which necessitates some special treatment. In short, it requires a special kind of grand canonical MC algorithm which I believe is worth sharing in this Note. Along the way, I reformulate the LGM from a different perspective, which may help better understand the physical meaning of the model.

II Materials and Method

II.1 Model

The basic structure of the LGM was presented in a previous paper Kinjo (2016). Here, I reformulate the model more formally with a few modifications. An LGM $\mathcal{M}$ is defined as a tuple of sets $\mathcal{M}=\{\mathcal{S},\mathcal{A},\mathcal{T}\}$. The set $\mathcal{S}$ is a set of $2N-1$ model sites, which is a disjoint union of the set of $N$ core sites $\mathcal{O}=\{O_{1},\cdots,O_{N}\}$ and the set of $N-1$ insert sites $\mathcal{I}=\{I_{1},\cdots,I_{N-1}\}$: $\mathcal{S}=\mathcal{O}\cup\mathcal{I}$. The number of core sites, $N$, is called the length of the model. The core sites represent those sequence positions that are present in a majority of sequences of the protein family of interest, and the insert sites represent the other positions. In this study, I define the core and insert sites as the positions of the “match” and “insert” states, respectively, of the profile hidden Markov model (HMM) of the corresponding Pfam Finn et al. (2014) family. The set $\mathcal{S}$ is partially ordered (denoted by “$\to$”) such that $O_{i}\to O_{i+1}$, $O_{i}\to I_{i}$, $I_{i}\to I_{i}$ and $I_{i}\to O_{i+1}$ for $i=1,\cdots,N-1$, and these ordered pairs are said to be potentially bonded (the dashed arrows in Fig. 1A). Those pairs of model sites that are not potentially bonded are called non-bonded. The set $\mathcal{A}$ is a set of sets of allowed amino acid residue types for each model sites: $\mathcal{A}=\{\mathcal{A}_{S}|S\in\mathcal{S}\}$ where $\mathcal{A}_{O_{i}}(i=1,\cdots,N)$ contains the 20 standard residues types and a symbol for the “delete” (“-”) and $\mathcal{A}_{I_{i}}(i=1,\cdots,N-1)$ contains the 20 standard residue types. The set $\mathcal{T}$ is a set of interacting non-bonded pairs of core sites. To define it precisely, we first define the (non-redundant) set of all pairs of core sites: ${\mathcal{T}^{*}}=\{(O_{i},O_{j})|O_{i},O_{j}\in\mathcal{O},i<j\}$. Then, $\mathcal{T}$ is defined as a subset: $\mathcal{T}\subset\mathcal{T}^{*}$. This subset is determined based on structural information (see below). A concrete example of a model with $N=6$ is shown in Fig. 1A.

Using this model, we can represent an alignment between an arbitrary amino acid sequence and the model. Let $\mathbf{a}=a_{1}\cdots a_{L}$ be an amino acid sequence of $L$ residues. An alignment between the sequence $\mathbf{a}$ and the model $\mathcal{M}$ is represented as a sequence of pairs of a model site (core or insert) and a residue: $\mathbf{X}=X_{1}\cdots X_{L_{\mathbf{X}}}$ where $L_{\mathbf{X}}$ is the length of the alignment and each $X_{i}$ is an aligned pair such as $(S,a)$ with $S\in\mathcal{S}$ and $a\in\mathcal{A}_{S}$, $a=a_{j}$ for some $j\in\{1,\cdots,L\}$ or $a=\text{``-''}$ (delete, only if $S$ is a core site). For $\mathbf{X}$ to be a proper alignment, there are two requirements. First, for any two consecutive aligned pairs $X_{k}X_{k+1}=(S,a)(S^{\prime},a^{\prime})$, the two model sites must be potentially bonded ($S\to S^{\prime}$). Note in particular that the same insert site can appear arbitrarily many times in an alignment due to the order $I_{i}\to I_{i}$. The two model sites that appear in two consecutive aligned pairs are called bonded (note the absence of the adverb “potentially”). Second, for any two aligned pairs $X_{k}=(S,a_{i})$ and $X_{l}=(S^{\prime},a_{j})$ in $\mathbf{X}$ that are aligned with proper amino acid residues $a_{i}$ and $a_{j}$ in the sequence $\mathbf{a}$, if $k<l$ then $i<j$. A concrete example of an alignment is shown in Fig. 1B.

In the present study, I have made a few modifications (Fig. 1) to the previous model Kinjo (2016). The first is that the alignment is now global with respect to the model but local with respect to the sequence (the previous model was global with respect to both the model and sequence). This means that the entire region of a model is always aligned with an amino acid sequence while only a part of the amino acid sequence may be aligned with the model (by ignoring possible flanking residues on both the N- and C-termini of the amino acid sequence).

Another major modification is the limited number of interacting non-bonded pairs (previously, all the core sites were interacting with each other: $\mathcal{T}=\mathcal{T}^{*}$). Between core sites more than 5 residues apart along the sequence, there may be interactions defined based on a representative native structure of the family. Two sites are defined to be interacting if the residues aligned to those sites are in contact in the corresponding (representative) native structure. Two residues are defined to be in contact if any non-hydrogen atoms in those residues are within 5Å. Interactions are defined only between core sites for simplicity.

Based on this representation of the alignment, the energy function of the alignment $\mathbf{X}$ is given as

[TABLE]

where $J$ and $K$ are, respectively, short-range and long-range interaction energy parameters to be determined from the given (observed) MSA and $\mu$ is the chemical potential, and $\mathcal{T}$ indicates the set of all interacting pairs of core sites and $s$ is a function to extract the model site from an aligned pair (i.e., $s(X)=S$ for $X=(S,a),S\in\mathcal{S}$). Thus, only $J$ and $K$ parameters constitute the intrinsic energy, and $\mu$’s are provided as external variables to control (perturb) the system.

We assume that the probability $P(\mathbf{X})$ of obtaining an alignment $\mathbf{X}$ is given by the Boltzmann distribution:

[TABLE]

where $T$ is the “temperature” and $\Xi[T]$ is the partition function

[TABLE]

Here, the summation is over all possible alignments ($\mathbf{X}$) with the model and all possible sequences. Since the alignment length can vary, this ensemble is considered to be a grand canonical ensemble. The grand potential is given by

[TABLE]

Given the parameters $J$, $K$ and $\mu$, we can sample sequences according to the probability distribution Eq. (2) by running Monte Carlo (MC) simulations.

The standard condition is defined to be the system with $\mu(X_{k})=0$ for all $X_{k}$, and the natural condition to be the standard condition with $T=1$.

The parameters $J$ and $K$ are determined iteratively so that the average residue pair counts for bonded and interacting non-bonded pairs over the samples produced by MC simulations under the natural condition match those observed in the given MSA of the family. However, this procedure is not straightforward since the length of alignment $\mathbf{X}$ is variable. The complication is due to the indexing scheme of aligned pairs where the meaning of each index, say “$k$” of $X_{k}$, is different for different alignments. Therefore, I reformulate the energy function by introducing some stochastic variables, single-site counts and bonded pair counts, based on model sites which are fixed for any alignments. In fact, this was the original formulation of the LGM in the previous paper Kinjo (2016).

For each model site and site pairs, we define the stochastic variables as functions of alignment $\mathbf{X}=X_{1}\cdots X_{L_{\mathbf{X}}}$. To do so, we first define the set $\mathcal{Y}$ of all possible pairs of model sites and amino acid residues as $\mathcal{Y}=\{(S,a)|S\in\mathcal{S},a\in\mathcal{A}_{S}\}$. Note that, while the set of all possible alignments $\{\mathbf{X}\}$ is an infinite set, the set $\mathcal{Y}$ is finite (namely, $|\mathcal{Y}|=21N+20(N-1)$). Now, the single-site count for the pair $Y\in\mathcal{Y}$ is defined as

[TABLE]

where $\delta_{Y,X_{k}}$ indicates Kronecker’s delta (i.e., $\delta_{Y,X_{k}}=1$ if $Y=X_{k}$ and $\delta_{Y,X_{k}}=0$ otherwise). Note that $n_{(O_{i},a)}(\mathbf{X})$ can be either 0 or 1 for any core site $O_{i}$ whereas $n_{(I_{i},a)}(\mathbf{X})$ for any insert site $I_{i}$ may have any values from 0 to infinity. Next, the bonded pair count for $Y,Y^{\prime}\in\mathcal{Y}$ is given as

[TABLE]

The bonded pair counts are defined only for potentially bonded pairs of model sites (those pairs of sites connected via dotted arrows in Fig. 1A). However, two neighboring model sites may not always be aligned with two consecutive residues in the amino acid sequence. For example, in the alignment $\mathbf{X}$ shown in Fig. 1B, we have the bonded pair count $n_{(O_{3},\mathtt{F}),(O_{4},\mathtt{-})}^{b}(\mathbf{X})=0$ because the pair of aligned sites $(O_{3},\mathtt{F})$ and $(O_{4},\mathtt{-})$ are not consecutive in the alignment $\mathbf{X}$ although the model sites $O_{3}$ and $O_{4}$ are potentially bonded in the model. Nevertheless, we have $n_{(O_{3},\mathtt{F})}(\mathbf{X})n_{(O_{4},\mathtt{-})}(\mathbf{X})=1$ because $n_{(O_{3},\mathtt{F})}(\mathbf{X})=1$ and $n_{(O_{4},\mathtt{-})}(\mathbf{X})=1$. Therefore, a bonded pair count cannot be reduced to a product of two single-site counts. In this manner, bonded pair counts account for the chain structure of polypeptide sequences, which should not be confused with the (quasi-)one-dimensional lattice structure of the model.

From the definitions of the single-site and bonded pair counts, we have the following relations Kinjo (2016). First, the single-site counts for each core site are normalized:

[TABLE]

Second, each single-site count is completely determined by bonded pair counts:

[TABLE]

Finally, it follows from the relation Eq. (7) that pair counts for non-bonded pairs also determine the single-site counts:

[TABLE]

These relations explain why chemical potentials $\mu$’s are not necessary as a part of the intrinsic energy parameters (Eq. 1). They also indicate that not all the variables are independent, which in turn indicates there is gauge freedom in the energy parameters $J$ and $K$. In this study, I (partially) fixed the gauge so that $J((O_{i},\mathtt{-}),(O_{i+1},\mathtt{-}))=0$ and $K((O_{i},\mathtt{-}),(O_{i+1},\mathtt{-}))=0$, that is, the interactions between the “delete” residues (“$\mathtt{-}$”) were defined to be zero.

Using the single-site counts and bonded pair counts, we can rewrite the energy function as

[TABLE]

where the summations $\sum_{Y,Y^{\prime}}^{b.p.}$ and $\sum_{Y,Y^{\prime}}^{n.b.p.}$ are taken over all the possible pairs of aligned pairs over potentially bonded pairs and interacting non-bonded pairs in the model sites, respectively.

II.2 Monte Carlo algorithm

I first describe the Metropolis-Hastings algorithm Iba (2005) for sampling alignments $\{\mathbf{X}\}$. Since alignment length $L_{\mathbf{X}}$ is variable depending on insertions, some care is necessary to ensure the detailed balance.

For a core site, we randomly select a residue $b$ out of 21 possible residue types, substitute it with the existing residue $a$ (Fig. 2A). Let the energies of the sequences before and after the substitution be $E_{a}$ and $E_{b}$, respectively, the substitution is accepted with the probability

[TABLE]

For an insert site, there are three possible moves: substituting a residue (Fig. 2B), extending the insert by adding a new residue or shortening the insertion by deleting a residue (Fig. 2C). Suppose an insert consists of $l$ residues ($l=0,1,\cdots$). Then there are $l+1$ possible positions for a new residue (one of 20 types is chosen with probability $q_{20}=1/20$) to be inserted between two consecutive residues, and $l$ possible positions (residues) for an existing residue to be deleted or substituted (deletion or a substituting residue is chosen with probability $q_{21}=1/21$). Accordingly, an extension of the insert is attempted with probability (attempt frequency) $p_{l}^{+}q_{20}$ where $p_{l}^{+}=(l+1)/(2l+1)$, and a deletion or substitution is attempted with probability $p_{l}^{-}q_{21}$ where $p_{l}^{-}=l/(2l+1)=1-p_{l}^{+}$. This choice of attempt frequencies automatically excludes the possibility of attempting a deletion or substitution when $l=0$.

For an extension of an insert, one of the $l+1$ positions and one of the 20 residue types are chosen randomly, then it is inserted to the selected position. Let the energies before and after the extension be $E_{l}$ and $E_{l+1}$, respectively. Then the extension is accepted with the probability

[TABLE]

where $p_{l+1}^{-}=(l+1)/(2l+3)$ (the probability of attempting substitution or deletion for an insert of length $l+1$), $q_{21}=1/21$ (uniform distribution for selecting a substituting residue or deletion) and $q_{20}=1/20$ (uniform distribution for selecting an inserted residue). The pre-exponential factor ensures the detailed balance condition Iba (2005).

For a substitution or deletion of an insert, randomly pick one of the $l$ positions in the insert. A deletion is attempted with probability $q_{21}=1/21$. Let the energy after the attempted deletion be $E_{l-1}$. Then the deletion is accepted with probability

[TABLE]

where $p_{l-1}^{+}=l/(2l-1)$ (the probability of attempting the extension of an insert of $l-1$ residues) and $p_{l}^{-}=1-p_{l}^{+}$. For a substitution, randomly pick a substituting residue with the uniform distribution ($q_{20}=1/20$), then apply the usual Metropolis criteria (c.f., Eq. 13).

For substitutions at core and insert sites, Gibbs sampling was also employed during parameter optimization. In these cases, substituting residues are selected according to the probability proportional to $\exp[-E(a)/T]$ where $E(a)$ is the energy of the sequence with the attempted substitution with residue type $a$ (either in a core or insert site).

For a model of length $N$, one sweep consists of $N$ and $N-1$ moves for randomly chosen core and insert sites, respectively.

II.3 Parameter optimization

Let $\braket{Q}_{\mathrm{obs}}$ denote the average value of the variable $Q(\mathbf{X})$ over a given MSA. If the MSA consists of $M$ aligned sequences ($\mathbf{X}^{1},\cdots,\mathbf{X}^{M}$), we compute the observed averages $\braket{n_{Y}}_{\mathrm{obs}}$, $\braket{n_{YY^{\prime}}^{b}}_{\mathrm{obs}}$ and $\braket{n_{Y}n_{Y^{\prime}}}_{\mathrm{obs}}$ as follow:

[TABLE]

where $q_{Y}$ (or $q_{Y^{\prime}}$) is 21 or 20 if the model site of the the aligned pair $Y$ (or $Y^{\prime}$) is a core site or an insert site, respectively, and $w_{t}$ is the position-based weight of the sequence $t$ Henikoff and Henikoff (1994). Only the core sites were used for computing the weights and the weights are normalized (i.e., $\sum_{t=1}^{M}w_{t}=1$). These average counts are referred to as number densities in the following. The value of the pseudocount $\gamma$ was set to a relatively small value ($\gamma=0.1$) in order to keep the expected sequence length with the pseudocounts closer to the value without the pseudocounts.

We define $\braket{Q}_{\mathrm{sim}}$ as the average of the variable $Q(\mathbf{X})$ over a set of samples obtained from MC simulations (at a constant temperature). It is given as a simple average:

[TABLE]

for a set of $M$ simulated samples.

Given the estimates for observed and simulated number densities, we can optimized the energy parameters $J$ and $K$ (Eqs. 1 and 12). Under the natural condition ($\mu(Y)=0$ for all $Y$ and $T=1$), the average energy of the observed MSA is given as

[TABLE]

The parameter optimization is done by maximizing $\braket{E}_{\mathrm{obs}}-\Omega[T]$ under the natural condition ($T=1$ in particular), or:

[TABLE]

where $\Omega[T]$ is the grand potential defined in Eq. (4). This may be regarded as a Legendre transform from the grand canonical ensemble with grand potential $\Omega[T]$ determined by $J(Y,Y^{\prime})$ and $K(Y,Y^{\prime})$ to the canonical ensemble with free energy $F[T]$ determined by the respective conjugate variables $\braket{n_{YY^{\prime}}^{b}}_{\mathrm{obs}}$ and $\braket{n_{Y}n_{Y^{\prime}}}_{\mathrm{obs}}$. The optimization of the parameters as given in this equation is equivalent to the principle of maximum entropy as used by others Morcos et al. (2011) and in the previous paper Kinjo (2016). Let $P$ denote either $J(Y,Y^{\prime})$ or $K(Y,Y^{\prime})$, and $Q(\mathbf{X})$ denote, respectively, $n_{YY^{\prime}}^{b}(\mathbf{X})$ or ${n_{Y}(\mathbf{X})n_{Y^{\prime}}(\mathbf{X})}$. We have the following relations:

[TABLE]

where $\braket{Q}$ is the exact average of $Q(\mathbf{X})$ obtained from the partition function, and approximated by the simulation average $\braket{Q}_{\mathrm{sim}}$. Based on these relations, we have the following procedure for parameter optimization.

Initialize $J$ and $K$ with some values. 2. 2.

Estimate $\braket{n_{YY^{\prime}}^{b}}_{\mathrm{sim}}$ and $\braket{n_{Y}n_{Y^{\prime}}}_{\mathrm{sim}}$ by running equilibrium MC simulations under the natural condition. 3. 3.

Update $P$ (either $J(Y,Y^{\prime})$ or $K(Y,Y^{\prime})$) by using the observed and simulated average values of $Q$ (respectively, $\braket{n_{YY^{\prime}}^{b}}$ or $\braket{n_{Y}n_{Y^{\prime}}}$)

[TABLE]

where $\alpha$ and $\beta$ are small parameters, and $R^{(\nu)}$ is a momentum term introduced to accelerate convergence Qian (1999). 4. 4.

Iterate 2 and 3.

The parameter $\alpha$ was set to 0.3, 0.1 or 0.01 depending on optimization stages (see Results), and $\beta$ was set to 0.95.

II.4 Data preparation

The multiple sequence alignment (MSA) and profile hidden Markov model (HMM) of the SH3 domain (Pfam PF00018) were downloaded from the Pfam database Finn et al. (2014) (version 30.0). The MSA was based on the “representative proteomes” of 75% sequence identity cutoff. Based on the profile HMM, the length of the model was set to 48.

A representative crystal structure was chosen for each family from the Protein Data Bank Berman et al. (2007) based on the criteria that there are no gaps (insert or delete) within the domain and the resolution is better than 2.0 Å. This is done by first querying the PDBj Mine2 relational database Kinjo et al. (2017), then by examining the alignments of the found PDB chains against the Pfam profile HMM using the HMMER hmmalign program Eddy (2011). Following this procedure, I selected chicken Src SH3 domain (PDB: 4HVU Bacarizo and Camara-Artigas (2013)) as a representative structure of the SH3 domain family, based on which there were 104 interacting non-bonded pairs of core sites.

III Results

III.1 Parameter optimization

I optimized the parameters $J$ and $K$ of the LGM model based on the MSA of representative sequences of the SH3 family (Pfam PF00018) under the natural condition ($T=1$ and $\mu(Y)=0$ for all $Y\in\mathcal{Y}$). The whole optimization process consisted of 3 stages. In the first stage, 120,000 sequences obtained from 12 trajectories of 10,000 sweeps each were used for computing the simulated densities, and then the parameters $J$ and $K$ were updated with $\alpha=0.3$ (Eq. 24). This optimization was repeated for 500 steps. In the second stage, 1,200,000 sequences from 12 trajectories of 100,000 sweeps each were for used for updating the parameters with $\alpha=0.1$. This stage consisted of 200 steps. The third stage was identical to the second stage except that $\alpha=0.01$.

By construction, a sufficiently optimized LGM reproduces the number densities of each model site as well as of bonded and (interacting) non-bonded pairs in the given MSA. This is indeed demonstrated in Fig. 3. For the single-site number densities, the correlation coefficient between observed and simulated values were greater than 0.99 ($>$0.9999 for the core sites, and 0.993 for the insert sites). The correlation is slightly weaker for the insert sites, suggesting the difficulty of sampling the arbitrarily large number of inserted residues (there can be infinite number of inserted residues in theory, but the number is always limited in simulations). Note that the single-site densities were not directly optimized, but they were optimized indirectly through the bonded and non-bonded pair densities (c.f., Eqs. 8-11). For bonded pairs, the correlation coefficients were $>$ 0.999, 0.993, 0.994 and 0.984 for $\braket{n_{(O_{i},a),(O_{i+1},b)}^{b}}$, $\braket{n_{(O_{i},a),(I_{i},b)}^{b}}$, $\braket{n_{(I_{i},a),(O_{i+1},b)}^{b}}$ and $\braket{n_{(I_{i},a),(I_{i},b)}^{b}}$, respectively. For the interacting non-bonded pair densities, the correlation was $>$0.999. We use these parameters in the following.

For comparison, I also examined the correlation between non-interacting non-bonded pairs which were left out from the optimization (Fig. 3D). The simulated values do correlate with the observed values, but, as expected, the correlation is not very high: 0.991 for $\braket{n_{O_{i}}(a)n_{O_{j}}(b)}$, $\sim$0.82 for both $\braket{n_{O_{i}}(a)n_{I_{j}}(b)}$ and $\braket{n_{I_{i}}(a)n_{I_{j}}(b)}$. The pairs involving insert sites are less correlated, indicating less sufficient sampling for insert sites.

III.2 Grand canonical MC simulations at various temperatures

Using the parameters defined above, MC simulations were performed at 3 different temperatures under the standard condition (i.e., $\mu(Y)=0$ for all aligned sites). At $T=1$ (the natural condition), the energy fluctuated mostly around 20-30 energy units (e.u.) with average 27.6 and standard deviation (s.d.) 6.4. The energy occasionally jumped to around 200 e.u., but very soon returned to the lower region (Fig. 4A, magenta line). At a high temperature ($T=1.1$), the energy steadily fluctuated around the average value of 239.2 e.u. with s.d. of 5.3 (Fig. 4A, green line). At an intermediate temperature of $T=1.02$, a two-state transition was observed (Fig. 4B) between low-energy and high-energy states.

One of the strengths of the LGM is that it can handle variable-length insertions. We define the sequence length of an alignment as the number of standard residues in the alignment (i.e., the total number of non-delete “residues”). The trajectories of sequence length show a large variety depending on the temperature (Figs. 4C,D). At $T=1$ and 1.1, the average sequence lengths were 48.8 (s.d. 4.6) and 91.5 (s.d. 6.9), respectively (Fig. 4C). The fluctuation of the sequence length is clearly correlated with that of the energy (Figs. 4A,C). A similar trend is also observed for $T=1.02$ (Fig. 4D).

To characterize the sequences generated at different temperatures, I performed homology searches using the hmmsearch program Eddy (2011) using the profile HMM of the SH3 domain (SH3_1.hmm provided by the Pfam database) against a database of sequences generated at a specified temperature. At $T=1.00$, 8,961 out of 10,000 (89.6%) of the sequences were significantly similar to the SH3 domain (E-value less than 0.01) whereas at $T=1.1$, no sequence was significantly similar. This confirms that not only the residue distribution on average, but also the MC-generated sequences at the natural conditions are similar to the natural sequences. At $T=1.02$, the fraction of significantly similar sequences was 41.3%. Thus, the transition observed in Fig. 4B indeed indicates an order-disorder transition in the sequence space. In the following, “ordered” or “disordered” sequences refer to those similar or dissimilar, respectively, to the natural sequences. Since the ordered sequences are likely to fold into the native fold of the SH3 domain whereas the disordered sequences are not, the transition may be regarded as a “folding” transition in the sequence space.

III.3 Perturbation on insert sites

In addition to changing the temperature, we can also perturb the system by introducing non-zero values for some or all of chemical potentials $\mu(Y)$. Here, I show the results of simulations with $\mu(I_{i},a)=\mu_{I}$, that is, the chemical potential was set to a constant value $\mu_{I}$ for all residue types $a$ for all insert sites $I_{i}$ with $T=1$ (Fig. 5). The case $\mu_{I}=0$ (with $T=1$) is the natural condition and the trajectory is identical to the one shown in Figs. 4A,C (magenta lines). When $\mu_{I}$ is set to a large value $\mu_{I}=0.2$, the energy (Fig. 5A) and the sequence length (Fig. 5C) mostly had large values as expected. At an intermediate value $\mu_{I}=0.1$, the trajectories again exhibited two-state transitions (Figs. 5B,D).

Out of the 10,000 sequences generated with $\mu_{I}=0.2$ and $\mu_{I}=0.1$, 2.5% and 46%, respectively, were significantly similar to the profile HMM of the SH3 domain. Since the chemical potentials at the core sites were kept to zero, this result indicates that prolonged insertions alone can trigger order-disorder transition.

III.4 Virtual temperature-jump experiment

Motivated by the existence of the “folding” transition, I next performed a virtual “temperature-jump” experiment in which the system was initially equilibrated at a high temperature ($T=1.1$) and then it was suddenly cooled to a lower temperature ($T=1.02$) and the trajectory was monitored (Fig. 6). As noted above (Fig. 4A), the system is disordered (non-natural-like) at $T=1.1$. An example, the alignment $\mathbf{X}_{270}$ (Fig. 6B), of the sequences shortly after the temperature jump illustrates some typical characteristics of high-energy, disordered sequences such as an excess of inserts and the lack of conserved residues (e.g., W32 and P47; c.f. $\mathbf{X}_{5885}$ below in Fig. 6B). As the “folding” proceeds, the sequence length tended to shrink ($\mathbf{X}_{1265}$ in Fig. 6B), and some conserved residues (e.g., W32) appeared when the energy became lower ($\mathbf{X}_{3964}$). However, the alignment still fluctuated to higher energies ($\mathbf{X}_{4153}$) before the entire region was ordered enough to significantly match the SH3 domain ($\mathbf{X}_{5885}$). Since $T=1.02$ is close to the transition temperature as suggested in the previous subsection (Fig. 4B), the alignment continued to fluctuate largely.

III.5 Mean-field approximation

One of the motivations for developing the MC method for the LGM model was that the mean-field approximation did not work well in the previous study Kinjo (2016) in which the parameters corresponding to $K$ in Eq. (1) were obtained by inverting the covariance matrix of the single-site counts of the core sites Morcos et al. (2011). This matrix inversion method has been derived from a mean-field approximation Plefka (1982); Kinjo (2015). Thus, the parameters based on the mean-field approximation was not consistent with the mean-field approximation of the LGM model. Note that the approximations are involved in two different situations, the first one in determining the parameters and the second one in computing the partition function. A detailed account for various approximations in the first case has been provided by Cocco et al. Cocco et al. (2017) Since the parameters are determined rigorously by MC simulations (within the sampling and numerical errors) in the present study, we can examine if it is indeed the case that the mean-field approximation of the LGM model is inappropriate in the second case. I applied the mean-field approximation to calculate the partition function of the LGM model using the parameters obtained from MC simulations. Please refer to the previous paper for the method to calculate the partition function using a transfer matrix method (with respective modifications to reflect the current model structure). Here, the mean-field of interactions is defined by

[TABLE]

where the summation is over the partners of the interacting non-bonded pairs. Note in particular that unlike the previous study Kinjo (2016), diagonal terms of the $K$ matrix are not included. The single-site number densities $\braket{n_{(O_{j},b)}}_{\mathrm{mf}}$ are given by

[TABLE]

where the probability of alignment $P_{\mathrm{mf}}(\mathbf{X})$ (c.f., Eq. 2) is computed using a transfer matrix method with the mean-field (Eq. 26) (see the previous paper Kinjo (2016) for the details). Eqs. (26) and (27) are mutually dependent, and therefore, are solved mutually consistently. For core sites, the difference of the mean-field or MC-generated single-site number densities from the observed ones was measured by the Kullback-Leibler (KL) divergence:

[TABLE]

where $\braket{n_{(O_{i},a)}}$ indicates either $\braket{n_{(O_{i},a)}}_{\mathrm{sim}}$ (obtained from MC simulations) or $\braket{n_{(O_{i},a)}}_{\mathrm{mf}}$ (obtained from the mean-field approximation).

The divergence of the core site number densities obtained from MC simulations were less than $10^{-3}$ for most of the sites (the only exception was $1.07\times 10^{-3}$ for the core site 29) whereas those obtained from the mean-field approximation were greater than 0.01 for 41 out of 48 core sites (and $>0.1$ for 15 sites) and as large as 0.82 for the site 46 (Fig. 7A). To check this was not an artifact of the exact partition function, I also optimized the $J$ parameters without long-range interactions so the truly exact partition function could be obtained (Fig. 7B). For this “$J$-only” system Kinjo (2016), the simulation and exact results were very similar and the divergence values were all small, mostly less than $2\times 10^{-3}$ (note that the “exact” result is not expected to yield divergence of 0 due to the errors in the parameters estimated from MC simulations).

For insert sites, we measure the difference between simulated and observed densities by the root square deviation instead of KL divergence (because the number densities are not normalized for insert sites):

[TABLE]

The deviation of the insert site number densities obtained from MC simulations were of order of $10^{-3}$, with the maximum value of $9\times 10^{-3}$ for $I_{40}$ whereas those obtained from the mean-field approximation were greater than 0.01 for 20 insert sites (Fig. 7C). For the “$J$-only” system, the deviations obtained from MC simulations and “exact” calculation were comparable (Fig. 7D). Since there are no non-bonded interactions involving insert sites, errors introduced to the core sites by the mean-field approximation also affect the insert sites through bonded interactions.

In summary, these results suggest that the mean-field approximation is indeed not appropriate for this model.

IV Discussion

In the present work, I used for long-range interactions only those involved in native 3D contacts. As such, the present model requires a priori knowledge of the native structure of the protein family of interest. While this may appear as a limitation in some respect (e.g., it cannot be used for structure prediction), it also has its advantages in other respects. First, parameter optimization is easier due to the smaller number of parameters. Since massive MC simulations are required to optimize a large number of parameters, this is clearly an advantage Lapedes et al. (1999). Second, discarding non-bonded pairs not involved in native contacts eliminates statistical noise irrelevant to the structural context. This enables a simpler interpretation of the long-range interactions. Moreover, it becomes also possible to investigate the role of the native structure in determining the conservation patterns of the family sequences by comparing with systems optimized with “misfolded” structures. This may be an interesting subject for a future study.

I have shown that the mean-field approximation gives inconsistent results for the LGM model. It is, of course, expected that any approximations yield different results than the exact solution to some extent, and the difference may or may not be acceptable depending on the application. For studying the patterns of sequence conservation, however, the degree of inconsistency resulting from the mean-field approximation (Fig.7A) is not acceptable. That is why in the previous study Kinjo (2016) I resorted to the Gaussian approximation which includes self interaction (diagonal) terms in addition to long-range interactions. The self interaction terms are intrinsic to each site and residue, which should not exist in principle except for possibly integrated-out effects of intermolecular interactions such as ligand binding. In other words, each residue at a certain position does not have a means to “know” where in the sequence it is located and how it should be conserved other than by interactions with other residues at other positions (or ligands). Thus, it seems more appropriate to employ MC simulations for the present purpose.

Using MC simulations, I have demonstrated the existence of a two-state transition in the sequence (sub)space around the SH3 domain. This sort of transition, analogous to the folding transition in the conformational space, has been suggested by Nishikawa Nishikawa (1993, 2002). If this type of transition is universal for many protein families, it has some interesting theoretical as well as practical implications. First, the boundary of each protein family can be determined clearly in terms of the transition point. Note that most of the conventional sequence models such as profile HMMs are essentially one-dimensional models where sharp transitions are simply impossible Landau and Lifshitz (1980). Therefore, the boundary between family members and non-members is necessarily fuzzy according to the conventional models. This limitation of the conventional methods may have already biased our knowledge of protein families. Second, if the boundary of a protein family in the sequence space can be clearly defined, it may be possible to assign artificially designed protein sequences to existing protein families by following mutational paths along which the ordered state is maintained. This means that proteins belong to the same family whether they are artificially designed or naturally selected. Third, it should be possible to characterize the transition state ensemble and to identify residues essential in determining protein families. Such characterization will be helpful in understanding evolutionary trajectories of proteins and what constitutes the protein family.

V Conclusion

In summary, Monte Carlo simulations of the LGM model will be a convenient means to explore the structure of the protein sequence space. More thorough investigations on the “folding” transition in the sequence space using this model is under way.

Acknowledgments

I thank Ikuo Fukuda and Ken Nishikawa for fruitful discussion.

Conflict of interest

None declared.

Author contributions

ARK did everything.