Biophysical models of cis-regulation as interpretable neural networks

The adoption of deep learning techniques in genomics has been hindered by the difficulty of mechanistically interpreting the models that these techniques produce. In recent years, a variety of post-hoc attribution methods have been proposed for addressing this neural network interpretability problem in the context of gene regulation. Here we describe a complementary way to address this problem. Our approach is based on the observation that two large classes of biophysical models for cis-regulatory mechanisms can be expressed as deep neural networks in which nodes and weights have explicit physiochemical interpretations. We also demonstrate how such biophysical networks can be rapidly learned, using modern deep learning frameworks, from the data produced by certain types of massively parallel reporter assays (MPRAs). These results suggest a scalable strategy for using MPRAs to systematically characterize the biophysical basis of gene regulation in a wide range of biological contexts. They also highlight gene regulation as an ideal venue for the development of scientifically interpretable approaches to deep learning.

Deep learning -the use of large multi-layer neural networks in machine learning applications -is revolutionizing information technology [1]. There is currently a great deal of interest in applying deep learning techniques to problems in genomics, especially for understanding gene regulation [2][3][4][5][6]. These applications of deep learning remain somewhat controversial, however, due to the difficulty of mechanistically interpreting trained neural networks.
Multiple strategies have been proposed for addressing this neural network interpretability problem [7][8][9][10][11][12][13]. A common feature of these attribution strategies is that they seek to extract meaning post-hoc from neural networks that have arbitrary architectures. However, there remains a substantial gap between the outputs of these attribution methods and fully mechanistic models of gene regulation.
Here we advocate for a complementary approach: the inference of neural network models whose architecture reflects explicit biophysical hypotheses for how a cis-regulatory sequence of interest might work. This strategy is based on two key observations. First, explicit biophysical models can be naturally formulated as deep neural networks in which nodes and weights have explicit physiochemical interpretations. This is true of standard thermodynamic models (which rely on a quasi-equilibrium assumption) [14][15][16][17][18][19] as well as fully kinetic models [20][21][22], and requires no mathematical approximations (c.f. [23,24]). Second, existing deep learning frameworks are able to rapidly infer such models from the data produced by certain classes of massively parallel reporter assays (MPRAs). Thermodynamic models are specified by a set of molecular complexes, or "states", which we index using s. Each state has both a Gibbs free energy ∆G s and an associated activity α s . These energies determine the probability P s of each state occurring in thermodynamic equilibrium via the Boltzmann distribution, 1 The energy of each state is, in turn, computed using integral combinations of the individual interaction energies ∆G j that occur in that state. We can therefore write ∆G s = j ω sj ∆G j , where ω sj is the number of times that interaction j occurs in state s. The resulting activity predicted by the model is given by the activities α s of the individual states averaged over this distribution, i.e., t = s α s P s . Fig. 1 illustrates a thermodynamic model for transcriptional activation at the E. coli lac promoter. This model involves two proteins, CRP and RNAP, as well as three interaction energies: ∆G C , ∆G R , and ∆G I . The rate of transcription t is further assumed to be proportional to the fraction of time that RNAP is bound to DNA (Fig. 1A). This model is summarized by four different states, two of which lead to transcription and two of which do not (Fig. 1B). Fig. 1C shows the resulting formula for t in terms of model parameters. This model is readily formulated as a feed-forward neural network (Fig. 1D). Indeed, all thermodynamic models of cis-regulation can be formulated as three-layer neural networks: layer 1 represents molecular interaction energies, layer 2 (which uses a softmin activation) represents state probabilities, and layer 3 (using linear activation) represents the biological activity of interest, which in this case is transcription rate.
We can infer thermodynamic models like these for a cis-regulatory sequence of interest (the wild-type sequence) from the data produced by a massively parallel reporter assay (MPRA) performed on an appropriate sequence library [25]. Indeed, a number of MPRAs have been performed with this explicit purpose in mind [25,[27][28][29][30]. To this end, such MPRAs are generally performed using libraries that consist of sequence variants that differ from the wild-type sequence by a small number of single nucleotide polymorphisms (SNPs). The key modeling assumption that motivates using libraries of this form is that the assayed sequence variants will form the same molecular complexes as the wild-type sequence, but with Gibbs free energies and state activities whose values vary from sequence to sequence. By contrast, variant libraries that contain insertions, deletions, or large regions of random DNA are unlikely to satisfy this modeling assumption. Fig. 2A summarizes the sort-seq MPRA described in [25]. Lac promoter variants were used to drive GFP expression in E. coli. GFP-expressing cells were then sorted into 10 bins using fluorescence-activated cell sorting, after which variant The parameter values inferred for the CRP energy matrix θ C , the RNAP energy matrix θ R , and the CRP-RNAP interaction energy ∆G I . Since increasingly negative energy corresponds to stronger binding, they y-axis in the logo plots is inverted. Logos were generated using Logomaker [26].
promoters within each bin were sequenced, yielding about 50K sequences in total. The authors then fit the biophysical model shown in Fig. 1C, but under the assumption that ∆G C = θ C · x + µ C and ∆G R = θ R · x + µ R , where x is a one-hot encoding of promoter DNA sequence.
For this study, we used TensorFlow to infer the same model formulated as a deep neural network. Specifically, we embellished the network in Fig. 1D by using same sequence-dependence for ∆G C and ∆G R as in [25]. To link t to the MPRA measurements, we introduced a feed-forward network with one hidden layer (having softmin activation) and a softmin output layer corresponding to the 10 bins into which cells were sorted. All model parameters were fit to the MPRA dataset using stochastic gradient descent and early stopping. The results agreed well with those reported in [25]. In particular, the parameters in the energy matrices for CRP ( θ C ) and RNAP ( θ R ) exhibited Pearson correlation coefficients of 0.986 and 0.994, respectively, with those reported in [25]. The protein-protein interaction energy that we found, ∆G I = −2.9 kcal/mol, was also compatible with the previously reported value ∆G I = −3.3 ± 0.4 kcal/mol.
A major difference between our results and those of [25] is the ease with which they were obtained. Training of the network in Fig. 2B consistently took about 15 minutes on a standard laptop computer. The model fitting procedure in [25], by contrast, relied on a custom Parallel Tempering Monte Carlo algorithm that took about a week to run on a multi-node computer cluster (personal communication), and more recent efforts to train biophysical models on MPRA data have encountered similar computational bottlenecks [29,30].
Also of note is the fact that in [25] the authors inferred models using information maximization. Specifically, the authors fit the parameters of t( x) by maximizing the mutual information between model predictions and observed bins I[t, bin].
One practical difficulty with this strategy is the need to estimate mutual information from finite data. Instead, we used maximum likelihood to infer the parameters of t( x) as well as the experimental transfer function (i.e., noise model) p(bin|t), which was modeled by a dense feed-forward network with one hidden layer. These two inference methods, however, are essentially equivalent: in the large data regime, the parameters of t that maximize I[t, bin] are the same as the parameters one obtains when maximizing likelihood over the parameters of both t and p(bin|t); see [31,32]. (B) A formula for the steady-state rate of mRNA production can be obtained using King-Altman diagrams [33,34]. (C) This formula is naturally represented using the three-layer neural network, where layer 1 represents transition rates, layer 2 represents the weights of distinct King-Altman diagrams, and layer 3 represents promoter activity. Here, black lines indicate weight 1, gray lines indicate weight 0, all activations are log-linear, and the dashed line represents L 1 normalization of layer 2 output.
A shortcoming of thermodynamic models is that they ignore non-equilibrium processes. Kinetic models address this problem, providing a fully non-equilibrium characterization of steady state activity. Kinetic models are specified by listing explicit state-to-state transition rates rather than Gibbs free energies. An example is shown in Fig. 3. King-Altman diagrams [33,34], a technique from mathematical enzymology, provide a straight-forward way to compute the steady-state occupancy of any individual state. Specifically, each state's occupancy is proportional to the sum of spanning trees that flow to that state, where each spanning tree's value is given by the product of rates comprising that tree. Fig. 3B illustrates this procedure for the kinetic model in Fig. 3A.
Here we have shown how thermodynamic and kinetic models of transcriptional regulation can be formulated as deep neural networks in which nodes and weights have explicit physiochemical meaning. We have further demonstrated how a thermodynamic model can be rapidly inferred from MPRA data using existing deep learning frameworks. These results suggest a new strategy for interpretable deep learning in the study of gene regulation, one complementary to existing post-hoc attribution methods.
Our approach can be applied to a wide variety of gene regulatory systems in both prokaryotes and eukaryotes. We demonstrated this approach in the context of a well-characterized bacterial promoter in order to avoid any uncertainty about the underlying biology, and because previous quantitative studies have established concrete results against which we could compare our inferred model. The same modeling approach, however, should be readily applicable to a wide variety of biological systems, including transcriptional regulation and alternative mRNA splicing in higher eukaryotes.
Some MPRA datasets are more amenable to this modeling strategy than others. Our approach is best-suited to mutagenesis studies of individual cis-regulatory sequences, rather than genome-wide MPRA datasets. By assaying SNP-containing variants of a specific cis-regulatory sequence, one preserves both the overall location of regulatory protein binding sites as well as the resulting protein-protein interactions. This greatly simplifies the biophysical models that one needs to infer. And indeed, such MPRA data have already been generated in a variety of eukaryotic systems [35,36,27,37]. The trade-off is that the resulting neural network model can only be expected to work in a narrow region of sequence space.
It may eventually be possible to relax this restriction while retaining biophysical interpretability via the use of convolutional or recurrent neural networks. Indeed, [24] has recently explored the possibility of using recurrent neural networks for biophysically modeling gene expression in flies. In principle, it should be possible to apply similar strategies to MPRAs performed on genome-wide or random sequence libraries [38,4,39,40,6]. But across all of biology, very few individual cis-regulatory sequences have been characterized to the level that our modeling strategy aims to elucidate. We therefore suggest that, at least in the short term, this approach should be used for modeling individual cis-regulatory sequences rather than genome-wide cis-regulatory codes.