Molecular Docking with Gaussian Boson Sampling

TL;DR

This paper demonstrates that Gaussian Boson Samplers can be used to predict molecular docking configurations, potentially aiding drug discovery by efficiently identifying stable binding poses in complex molecular systems.

Contribution

It introduces a novel approach using Gaussian Boson Samplers to solve the molecular docking problem by sampling large-weight cliques in a graph, improving prediction of binding configurations.

Findings

Gaussian Boson Samplers can sample stable docking configurations with high probability.

The approach enhances classical algorithms for molecular docking.

Benchmarking on a target enzyme shows promising results.

Abstract

Gaussian Boson Samplers are photonic quantum devices with the potential to perform tasks that are intractable for classical systems. As with other near-term quantum technologies, an outstanding challenge is to identify specific problems of practical interest where these quantum devices can prove useful. Here we show that Gaussian Boson Samplers can be used to predict molecular docking configurations: the spatial orientations that molecules assume when they bind to larger proteins. Molecular docking is a central problem for pharmaceutical drug design, where docking configurations must be predicted for large numbers of candidate molecules. We develop a vertex-weighted binding interaction graph approach, where the molecular docking problem is reduced to finding the maximum weighted clique in a graph. We show that Gaussian Boson Samplers can be programmed to sample large-weight cliques, i.e., stable docking configurations, with high probability, even in the presence of photon loss. We also describe how outputs from the device can be used to enhance the performance of classical algorithms and increase their success rate of finding the molecular binding pose. To benchmark our approach, we predict the binding mode of a small molecule ligand to the tumor necrosis factor-${\alpha}$ converting enzyme, a target linked to immune system diseases and cancer.

Figures

Tables

Table 1. Table S1: Knowledge-based pharmacophore potential . Data is derived from the PDBbind dataset from 2015 Liu  et al. ( 2014 ); Wang  et al. ( 2005 , 2004 ) . The matrix is lower-diagonal since any potential function is symmetric.

Pharmacophore type	Negative charge	Positive charge	Hydrogen-bond donor	Hydrogen-bond acceptor	Hydrophobe	Aromatic
Negative charge	0.2953
Positive charge	0.6459	0.1596
Hydrogen-bond donor	0.7114	0.4781	0.5244
Hydrogen-bond acceptor	0.6450	0.7029	0.6686	0.5478
Hydrophobe	0.1802	0.0679	0.1453	0.2317	0.0504
Aromatic	0.0	0.1555	0.1091	0.0770	0.0795	0.1943

Equations

$$p(S)=\frac{\mathrm{Haf}({\mathcal{A}}_{\mathrm{S}})}{n_1!\dots n_M!\sqrt{\det (\sigma +\text{mathbbm}1/2)}},$$

$$\mathrm{Haf}(A)=\sum _{\mathcal{M}\in \mathrm{PMP}}\prod _{(i,j)\in \mathcal{M}}{A}_{ij},$$

$$B=\Omega (D-A)\Omega ,$$

$$p(S)\propto [\det ({\Omega }_S)\mathrm{Haf}(A_S){]}^2,$$

$${\Omega }_{ii}=c(1+\alpha w_i),$$

$$\underset{f}{\sup }\frac{f^{T}Bf}{f^{T}f}$$

$$=\underset{g}{\sup }\frac{g^T\tilde{L}g}{g^{T}g}\frac{g^{T}g}{g^T{D}^{-1/2}{\Omega }^{-2}{D}^{-1/2}g}\le$$

$$\le \underset{g}{\sup }\frac{2g^{T}g}{g^T{D}^{-1/2}{\Omega }^{-2}{D}^{-1/2}g}=$$

$$=\underset{f}{\sup }\frac{2f^T\Omega D\Omega f}{f^{T}f}\le \underset{j}{\max }2c^2{\omega }_j^2{d}_j,$$

$$\mathrm{Haf}(B)=\det (\Omega )\mathrm{Haf}(A).$$

$$\mathrm{Haf}(B)$$

$$\stackrel{(b)}{=}\sum _{M\in \mathrm{PMP}}\prod _{(ij)\in M}{A}_{ij}\prod _{(ij)\in M}({\omega }_i{\omega }_j)=$$

$$\stackrel{(c)}{=}\prod _k({\omega }_k)\sum _{M\in \mathrm{PMP}}\prod _{(ij)\in M}{A}_{ij},$$

$$p(S_N)=\mathrm{Tr}\left[\prod _{i\in S_N}{P}_1^i\prod _{i\notin S_N}{P}_0^i\rho \right],$$

$$N=\sum _{j=1}^M\mathrm{Tr}[P_1^j\rho ]=M-\sum _{j=1}^M⟨0|{\rho }_j|0⟩,$$

$$N[\sigma ]=M-\sum _{j=1}^M\frac{1}{\sqrt{\det [{\sigma }_j+\text{openone}/2]}},$$

$$V_{S}=\big{\{}v_{i}\,|\,\bm{p_{i}}\in S\big{\}}$$

$$E_{S}=\big{\{}(v_{i},v_{j})\,|\,\bm{p_{i},p_{j}}\in S,i<j\big{\}}$$

$${\omega }_S((v_i,v_j))=||{\mathbf{p}}_{\mathbf{i}}-{\mathbf{p}}_{\mathbf{j}}||$$

$${\alpha }_S:V_S\to \mathbb{L}$$

$$I_{L,B}=(\mathcal{V},\mathcal{E},\kappa ,\tau ,ϵ),$$

$$\mathcal{V}=V_L\times V_B,$$

$$\mathcal{E}\subseteq\Big{\{}\big{(}(v_{l_{1}},v_{b_{1}}),(v_{l_{2}},v_{b_{2}})\big{)}\,|\,v_{l_{1}},v_{l_{2}}\in V_{L},v_{b_{1}},v_{b_{2}}\in V_{B}\big{\}}$$

$$\Omega\big{(}(v_{l},v_{b})\big{)}=\kappa\big{(}\alpha_{L}(v_{l}),\alpha_{B}(v_{b})\big{)},$$

$$P_{\text{refl}}(i,j)=\max (P)-\min (P)-P_{\text{orig}}(i,j),$$

Full text

I Introduction

In his lecture “Simulating Physics with Computers” Feynman (1982), Richard Feynman famously argued that classical computing techniques alone are insufficient to simulate quantum physics. Since then, significant progress has been made in formalizing this intuition by finding explicit examples of quantum systems whose classical simulation can be convincingly shown to require exponential resources. An example is Boson Sampling, first introduced by Aaronson and Arkhipov Aaronson and Arkhipov (2011). In this paradigm, identical photons interfere by passing through a network of beam-splitters and phase-shifters, and are subsequently detected at the output ports of the network. Despite the simplicity of this model, it has been shown that, under standard complexity-theoretic conjectures, generating samples from the output photon distribution requires exponential time on a classical computer Aaronson and Arkhipov (2011); Clifford and Clifford (2018); Neville et al. (2017). Several variants of boson sampling have been proposed that aim at decreasing the technical challenges with its experimental implementation Lund et al. (2014); Bentivegna et al. (2015); Latmiral et al. (2016); Hamilton et al. (2017); Chabaud et al. (2017); Lund et al. (2017); Hamilton et al. (2017); Kruse et al. (2018); Vernon et al. (2018).

Most efforts in the study of Boson Sampling have been focused on its viability to disprove the Extended Church-Turing thesis Wolfram (1985); not on its potential practical applications. Nevertheless, it is possible to ask: if Boson Sampling devices are powerful enough that they cannot be simulated with conventional computers, is there a way of programming them to perform a useful task? In fact, practical applications of Boson Sampling have already been reported. In Ref. Huh et al. (2015), it was shown that a Boson Sampling device can be used to efficiently estimate the vibronic spectra of molecules, a problem for which in general no efficient algorithm is known. Proof-of-principle demonstrations have also been reported Clements et al. (2017); Sparrow et al. (2018). Additionally, Refs. Arrazola and Bromley (2018); Arrazola et al. (2018); Brádler et al. (2018) discuss how a specific model known as Gaussian Boson Sampling (GBS) can be employed in combinatorial optimization problems concerned with identifying large clusters of data.

Molecular docking is a computational method for predicting the optimal interaction of two molecules, typically a small molecule ligand and a target receptor. This method works by searching the configurational space of the two molecules and scoring each pose using a potential energy function. Using molecular structures to determine stable ligand-receptor complexes is a central problem in pharmaceutical drug design De Vivo et al. (2016); Kitchen et al. (2004); Meng et al. (2011); Gschwend et al. (1996); Shoichet et al. (2002). Several techniques for finding stable ligand-receptor configurations have been developed, including shape-complementarity methods Shentu et al. (2008); Shoichet et al. (1992); Shoichet and Kuntz (1993); Goldman and Wipke (2000); DesJarlais et al. (1986); Dias et al. (2008) and molecular simulation of the ligand-receptor interactions Wu et al. (2003); Alonso et al. (2006), which vary in their computational requirements. For high-throughput virtual screening of large chemical libraries, it is desirable to search and score ligand-receptor configurations using as few computational resources as possible Shoichet (2004). Motivated by these computational problems, several recent efforts have focused on practical applications of near-term quantum computers in the life sciences Moll et al. (2017); Streif et al. (2018); O’malley et al. (2016); Reiher et al. (2017); Fingerhuth et al. (2018); Babej et al. (2018); Li et al. (2018); Hernandez and Aramon (2017); Hernandez et al. (2016).

In this work, we show how GBS can be used to solve the molecular docking problem. We extend the binding interaction graph approach, where the problem of identifying docking configurations can be reduced to finding large clusters in weighted graphs Kuhl et al. (1984); Kuntz et al. (1982). We then show how GBS devices can be programmed to sample from distributions that assign large probabilities to these clusters, thus helping in their identification. Docking configurations can be obtained by direct sampling or by hybrid algorithms where the GBS outputs are post-processed using classical techniques. We apply our method through numerical simulations to find molecular docking configurations for a known ligand-receptor interaction Rao et al. (2007). Several therapeutic agents targeting this protein have entered into clinical trials for both cancer and inflammatory diseases Moss et al. (2008a).

II Background

Before presenting our results we provide relevant background information on Gaussian Boson Sampling, graph theory, and molecular docking.

II.1 Gaussian Boson Sampling

Quantum systems such as the quantum harmonic oscillator or the quantized electromagnetic field can be described by phase-space methods. Here, each state is uniquely determined by a quasi-probability distribution such as the Wigner function $W(x,p)$ over its position $x$ and momentum $p$ variables Serafini (2017). A quantum state is called Gaussian if its Wigner function is Gaussian Weedbrook et al. (2012). Any multi-mode Gaussian state $\rho$ is parametrized by its first and second moments, namely the displacement $\alpha_{j}=\operatorname{Tr}[\rho\hat{\xi}_{j}]$ and the covariance matrix $\sigma$ with entries $\sigma_{jk}=\operatorname{Tr}[\{\hat{\xi}_{j},\hat{\xi}_{k}\}]/2$, where $\hat{\xi}_{j}$ is a vector of creation and annihilation operators: calling $M$ the number of modes, $\hat{\xi}_{j}=\hat{a}_{j}=(\hat{x}_{j}+i\hat{p}_{j})/\sqrt{2}$ and $\hat{\xi}_{M+j}=\hat{a}_{j}^{\dagger}$ for $j=1,\dots,M$. Gaussian quantum states are ubiquitous in quantum optics, and have enabled detailed theoretical modeling and coherent manipulations in experiments Weedbrook et al. (2012); Vogel and Welsch (2006).

In spite of their infinite-dimensional Hilbert space, Gaussian states can be simulated efficiently, as their evolution can be modeled by linear transformations such as Bogoliubov rotations Bartlett et al. (2002). However, when non-Gaussian measurements are employed, e.g., via photon-counting detectors Lund et al. (2014); Hamilton et al. (2017) or threshold detectors Quesada et al. (2018), modelling measurement outcomes becomes extremely challenging even for supercomputers. Indeed, it has been shown that under standard complexity assumptions, sampling from the resulting probability distribution cannot be done in polynomial time using classical resources Aaronson and Arkhipov (2011); Hamilton et al. (2017); Kruse et al. (2018).

For a Gaussian state with zero displacement and covariance matrix $\sigma$, the Gaussian Boson Sampling (GBS) distribution obtained by measuring the state with photon-counting detectors is given by Hamilton et al. (2017):

[TABLE]

where $\mathcal{A}=\begin{pmatrix}0&\operatorname{\mathbbm{1}}\\ \operatorname{\mathbbm{1}}&0\end{pmatrix}\left[\operatorname{\mathbbm{1}}-(\sigma+\operatorname{\mathbbm{1}}/2)^{-1}\right]~{},$ and $\mathcal{A}_{S}$ is a $2N\times 2N$ submatrix of $\mathcal{A}$, with $N=\sum_{j=1}^{M}n_{j}$. The set $S=(n_{1},\dots,n_{M})$ defines a measurement outcome, where $n_{j}$ is the number of photons in mode $j$, and the submatrix $\mathcal{A}_{S}$ is obtained by selecting rows and columns of $\mathcal{A}$, as described in Ref. Hamilton et al. (2017). The function ${\rm Haf}(\mathcal{A}_{S})$ is the Hafnian of $\mathcal{A}_{S}$, a matrix function which is #P-Hard to approximate for worst-case instances Caianiello (1953); Barvinok (2016); Björklund et al. (2018). For a $2N\times 2N$ matrix $A$, it is defined as

[TABLE]

where $\rm PMP$ is the set of perfect matching permutations, namely the possible ways of partitioning the set ${1,\dots,2N}$ into subsets of size 2. When threshold detectors are employed Quesada et al. (2018), the output is a binary variable $s_{j}$ for each mode: $s_{j}=1$ corresponds to a “click” from the $j$th detector that occurs whenever $n_{j}>0$; on the other hand, $s_{j}=0$ for $n_{j}=0$. The probability distribution with threshold detectors can be obtained by summing infinitely many probabilities from Eq. (1) or via closed-form expressions that require evaluating an exponential number of matrix determinants Quesada et al. (2018).

II.2 GBS to find dense subgraphs

When $A$ is the adjacency matrix of an unweighted graph $G$, the Hafnian of $A$ is equal to the number of perfect matchings in $G$. Using mathematical properties of the Hafnian, it was shown in Ref. Brádler et al. (2018) that a GBS device can be programmed to sample from a distribution $p(S)\propto\frac{|{\rm Haf}(A_{S})|^{2}}{c^{N}}$. The parameter $c$ depends on the spectral properties of $A$ and can be tuned to lower the probability of observing photon collisions, i.e., $n_{j}\geq 2$ for some $j$. More details are provided in Appendix A. In the collision-free subspace, $A_{S}$ is the adjacency matrix of the subgraph specified by the vertices $j$ for which $n_{j}=1$, and ${\rm Haf}(A_{S})$ is equal to the number of perfect matchings in this subgraph. Therefore, a GBS device can be programmed to sample large-Hafnian subgraphs with high probability.

The density of a graph $G$ is defined as the number of edges in $G$ divided by the number of edges of the complete graph. Intuitively, a subgraph with a high number of perfect matchings should have a large density; a connection that was made rigorous in Ref. Aaghabali et al. (2015). This fact was used in Ref. Arrazola and Bromley (2018) to show that GBS devices can be programmed to sample dense subgraphs with high probability. Hybrid quantum-classical optimization algorithms can be built by combining GBS random sampling with stochastic optimization algorithms for dense subgraph identification.

II.3 Molecular docking

Molecular docking is a computational tool for rational structure-based drug discovery. Docking algorithms predict non-covalent interactions between a drug molecule (ligand) and a target macromolecule (receptor) starting from unbound three-dimensional structures of both components. The output of such algorithms are predicted three-dimensional orientations of the ligand with respect to the receptor binding site and the respective score for each orientation. Reliable determination of the most probable ligand orientation, and its ranking within a series of compounds, requires accurate scoring functions and efficient search algorithms Halperin et al. (2002). The scoring function contains a collection of physical or empirical parameters that are sufficient to score binding orientation and interactions in agreement with experimentally determined data on active and inactive ligands. The search algorithm describes an optimization approach that can be used to obtain the minimum of a scoring function, typically by scanning across translational and rotational degrees of freedom of the ligand in the chemical environment of the receptor. In the simplest case, both the ligand and the receptor can be approximated as rigid bodies, but more accurate methods can account for inherent flexibility of the ligand and receptor Dias et al. (2008). As is the case for most molecular modelling approaches, a trade-off exists between accuracy and speed.

High-performance algorithms enable molecular docking to be used for screening large compound libraries against one or more protein targets. Molecular docking and structure-based virtual screening are routinely used in pharmaceutical research and development Kalyaanamoorthy and Chen (2011). However, evaluating billions of compounds requires accurate and computationally efficient algorithms for binding pose prediction. Widely used approaches for molecular docking employ heuristic search methods (simulated annealing Yue (1990) and evolutionary algorithms Trott and Olson (2009)) and deterministic methods Rarey et al. (1996). In one combinatorial formulation of the binding problem utilized in the DOCK 4.0 and FLOG Ewing et al. (2001); Miller et al. (1994), an isomorphous subgraph matching method is utilized to generate ligand orientations in the binding site Kuntz et al. (1982); Kuhl et al. (1984); Ewing and Kuntz (1997). In this formulation of the binding problem, both the ligand and the binding site of the receptor are represented as complete graphs. The vertices of these graphs are points that define molecular geometry and edges capture the Euclidean distance between these points. In order to strike a balance between the expressiveness of the graph and its size, we reduce the all-atom molecular models of the ligand and receptor to a pharmacophore representation Güner (2000); Yang (2010).

A pharmacophore is a set of points which have a large influence on the molecule’s pharmacological and biological interactions. These points may define a common subset of features, such as charged chemical groups or hydrophobic regions, that may be shared across a larger group of active compounds. For the purposes of this study, we define six different types of pharmacophore points: negative/positive charge, hydrogen-bond donor/acceptor, hydrophobe, and aromatic ring. In the graph representation, the type of the pharmacophore point is preserved as a label associated with its vertex. Hence we refer to this molecular graph representation as a labeled distance graph (see also Appendix B). As illustrated in Fig. 1, a labeled distance graph is constructed as follows for both the ligand and receptor:

Heuristically identify pharmacophore points likely to be involved in the binding interaction. These form the vertices of the graph. 2. 2.

Add an edge between every pair of vertices and set its weight to the Euclidean distance between the pharmacophore points they represent. 3. 3.

Assign a label to every vertex according to the respective type of pharmacophore point it represents.

III GBS for molecular docking

III.1 Mapping molecular docking to maximum weighted clique

The labeled distance graphs described in Section II.3 capture the geometric three-dimensional shapes and the molecular features of both the protein binding site and the ligand that interacts with it. In this section, akin to Kuhl et al. (1984), we combine these two graphs into a single binding interaction graph. Subsequently, we reduce the molecular docking problem to the problem of finding the maximum weighted clique.

If two pharmacophore points are interacting, they form a contact. A binding pose can be defined by a set of three or more contacts that are not colinear. We model contacts as pairs of interacting vertices of the labeled distance graphs of the ligand and the binding site. Consider the labeled distance graph $G_{L}$ of the ligand and the labeled distance graph $G_{B}$ of the binding site, with their vertex sets $V_{L}$ and $V_{B}$ respectively. A contact is then represented by a singe vertex $c_{i}\in V_{L}\times V_{B}$. The set of possible contacts forms the vertices of the binding interaction graph. In principle, any pharmacophore point of the ligand could be interacting with any pharmacophore point of the binding site, and therefore we have to consider every possible pair of corresponding interacting vertices. Hence the number of vertices of the binding interaction graph is $nm$, where $n$ is the number of vertices of the labeled distance graph $G_{L}$ and $m$ is the number of vertices of the labeled distance graph $G_{B}$.

The goal of the binding interaction graph is to model possible binding poses via sets of contacts. However, not every combination of contacts is physically realizable. Two contacts are not be compatible if their mutual realization would violate the geometrical shapes of the ligand and the binding site. To model this, the binding interaction graph contains an edge between two contacts if and only if they are compatible. As a result, a pairwise compatible set of contacts, i.e., such as would arise from a true binding pose, forms a complete subgraph of the binding interaction graph. A complete subgraph, also called a clique, in a graph $G$ is a subgraph where all possible pairs of vertices are connected by an edge.

The compatibility of contacts is captured by the notion of $\tau$ flexibility, which is illustrated in Fig. 2 (see also Appendix B). Even though both the ligand and the binding site can exhibit a certain amount of flexibility, in general, geometric distances between two contacts have to be approximately the same both on the ligand and the binding site. Two contacts $(v_{l_{1}},v_{b_{1}})$ and $(v_{l_{2}},v_{b_{2}})$ form a $\tau$* flexible contact pair* if the distance between the pharmacophore points on the ligand (points corresponding to vertices $v_{l_{1}}$ and $v_{l_{2}}$) and the distance between the pharmacophore points on the binding site (points corresponding to vertices $v_{b_{1}}$ and $v_{b_{2}}$) does not differ by more than $\tau+2\epsilon$ (see Panel C in Fig. 2). The constants $\tau$ and $\epsilon$ describe the flexibility constant and interaction distance respectively.

In order to model varying interaction strengths between different types of pharmacophore points, we associate a different weight to every vertex of the binding interaction graph. The weights are derived using the pharmacophore labels that are captured in the labeled distance graphs of the ligand and the binding site. Given a set of labels $\mathbb{L}$, a potential function $\kappa:\mathbb{L}\times\mathbb{L}\rightarrow\mathbb{R}$ is applied to compute the weights of the individual vertices. This allows us to bias the algorithm towards stronger intermolecular interactions. Potential functions can be derived in several ways, ranging from pure data-based approaches such as statistical or knowledge-based potentials Poole and Ranganathan (2006); Mintseris et al. (2007); Gohlke and Klebe (2001) to quantum-mechanical potentials Kitchen et al. (2004). Details of the potential used in this study are described in Section IV.

Hence under the model derived in this study, the most likely binding poses correspond to vertex-heaviest cliques in the binding interaction graph. The problem of finding a maximum weighted clique is a generalization of the maximum clique problem of finding the clique with the maximum number of vertices. When $G$ has $n$ vertices, the number of possible subgraphs is $\mathcal{O}(2^{n})$, so a brute force approach becomes rapidly infeasible for growing values of $n$. The max-clique decision problem is NP-hard Karp (1972): as such, unless P=NP, in the worst case any exact algorithm run for superpolynomial time before finding the solution. There are deterministic and stochastic classical algorithms for finding both the maximum cliques and maximum weighted cliques, or for finding good approximations when $n$ is large Wu and Hao (2015).

III.2 Max weighted clique from GBS

In this section, we show that a GBS device can be programmed to sample from a distribution that outputs the max-weighted clique with high probability. The main technical challenge is to program a GBS device to sample, with high probability, subgraphs with a large total weight that are as close as possible to a clique. Consider the graph Laplacian $L=D-A$, where $D$ is the degree matrix and $A$ the adjacency matrix. The normalized Laplacian Chung and Graham (1997) $\tilde{L}=D^{-1/2}LD^{-1/2}$ is positive semidefinite and its spectrum is contained in $[0,2]$. More generally, we define a rescaled matrix

[TABLE]

where $\Omega$ is a suitable diagonal matrix. If the largest entry of $\Omega$ is bounded as shown in Appendix A, then the spectrum of $B$ is contained in $[0,c]$, where $c\leq 1$ can be tuned depending on the maximum amount of squeezing obtainable experimentally. Using the decoupling theorem from Appendix A, we find that a GBS device can be programmed to sample from the distribution

[TABLE]

where we consider outputs $S=(n_{1},\ldots,n_{M})$ with $n_{j}\leq 1$ and $N=\sum_{j}n_{j}$ total photons. In the collision-free subspace, the dependence on the diagonal matrix $D$ disappears so we may focus on programming GBS with a rescaled adjacency matrix $\Omega A\Omega$. From a GBS sample $S$, we construct the subgraph $H$ of $G$ made by vertices $j$ with $n_{j}=1$. The matrix $A_{S}$ is the $N\times N$ adjacency matrix of $H$. The Hafnian of an adjacency matrix is maximum for the complete graph, namely when $H$ is a clique. Therefore, for a fixed total number of photons $N$, the Hafnian term maximizes the probability of detecting photon configurations that correspond to a clique.

Different choices are possible for the weighting matrix $\Omega$. For an unweighted graph, convenient choices are either a constant $\Omega$ or $\Omega\propto D$. In the former case, $\det\Omega_{S}=c^{N}$ for $c<1$, so the parameter $c$ can be tuned via squeezing in order to penalize larger $N$, i.e., larger subgraphs (see Appendix A.3). In the latter case, $\det\Omega=c^{N}\det D$ is proportional to the Narumi-Katayama index Gutman and Ghorbani (2012), which describes some topological properties of the graph. Similarly to the Hafnian, it is maximum when $H$ is a clique.

For a vertex-weighted graph, we can use the freedom of choosing $\Omega$ to favour subgraphs with larger total weight. There are multiple ways of introducing the weights $w_{j}$ in $\Omega$ and a convenient choice is

[TABLE]

where $c$ is a normalization to ensure the correct spectral properties and $\alpha>0$ is a constant. When $\alpha$ is small, the determinant term $\det\Omega_{S}\approx 1+\alpha\sum_{j:n_{j}=1}w_{j}$ is large when the subgraph $H$ has a large total weight. This is useful for the max-weighted clique problem as it introduces a useful bias in the GBS probability of Eq. (4) that favours heavier subgraphs. However, if $\alpha$ is too large, the Hafnian term in Eq. (4) becomes less important and GBS will sample heavy subgraphs that typically do not contain cliques. To prevent this occurrence, the parameter $\alpha$ must be chosen carefully. Ideally, the weights should give a positive bias to heavy cliques, but should not favour heavy subgraphs that are not cliques. More details are discussed in Appendix A.

III.3 Hybrid algorithms

GBS devices can in principle have a very high sampling rate – primarily limited by detector dead time – so just by observing the photon distribution it is possible to extract the maximum weighted clique for small enough graphs. We call this simple strategy GBS random search – see Fig. 3 for a graphical explanation of the method. However, selecting photon outcomes that correspond only to cliques means wasting samples that are potentially close to the solution. Indeed, an optimally programmed GBS device will sample from both the correct solution and neighboring configurations with high probability. Therefore, we propose two algorithms to post-process all GBS data which incur an overhead in run time but are especially useful for finding cliques in larger graphs.

Greedy Shrinking: Starting from an output subgraph $H$ from GBS, vertices are removed based on a local rule until a clique is found – see Fig. 3 for a graphical explanation of the method. Removal is based on vertex degree and weight. Vertices with small degree are unlikely to be part of a clique making them good candidates to be discarded. The role of the weights is less straightforward: vertices with low weight may not be part of the max-weighted clique, but this assumption may be incorrect if the clique is made by a heavy core together with a few light vertices. Because of this, vertex degree is prioritized over vertex weight during the greedy shrinking stage. More precisely, the algorithm proceeds as follows:

From a GBS outcome, build a subgraph $H$ with vertices corresponding to the detectors that “click”. 2. 2.

If $H$ is a clique, return $H$. 3. 3.

Otherwise, set $v$ as the set of vertices in $H$ with smallest degree. 4. 4.

Set $w$ as the subset of $v$ with lowest weight. 5. 5.

Remove a random element of $w$ from $H$ and go back to step 2.

Expansion with Local Search: GBS provides high-rate samples from max-cliques, and greedy shrinking enhances the probability of finding a solution via classical post-processing of sampled configurations. We may increase the probability of finding the solution even further, at the cost of a few more classical steps. This is done by employing a local search algorithm that tries to expand the clique with neighbouring vertices, as shown also in Fig. 3. Algorithms such as Dynamic Local Search (DLS) Pullan and Hoos (2006) and Phased Local Search (PLS) Pullan (2006) are among the best-performing classical algorithms for max-clique Wu and Hao (2015). These algorithms usually start with a candidate clique formed by a single random vertex, and then try to expand the clique size and replace some of its vertices by locally exploring the neighbourhood. More precisely, the following iteration is repeated until a sufficiently good solution is found, or the maximal number of steps is reached:

Grow stage: Starting from a given clique, generate the set of vertices that are connected to all vertices in the clique. If this set is non-empty, select one vertex at random, possibly with large weight, and add it to the clique. 2. 2.

Swap stage: If the above set is empty, generate the set of vertices that are connected to all vertices in the clique except one (say $v$). From this new set, select a vertex at random and swap it with $v$. This gives a new clique of the same size but with different vertices, thus constituting a local change to the clique. For max-weighted clique, the swapping rule also considers vertex weight.

An important aspect of the above local search is that, at each iteration step, the candidate solution is always a clique and the algorithm tries to expand it as much as possible. GBS can be included in this strategy in view of its ability to provide a starting configuration that is not a mere random vertex. Indeed, a GBS output after greedy shrinking is always a clique, with a comparatively large probability of being close to the maximum clique. In case the candidate output from greedy shrinking is not the maximum clique, then it can be expanded with a few iterations of local search. Since the cliques sampled from a carefully programmed GBS device are, with high probability, larger than just a random vertex, the number of classical expansion steps is expected to be significantly reduced. This will be demonstrated with relevant numerical examples in the following section.

IV Numerical results

We study the binding interaction between the tumor necrosis factor-$\alpha$ converting enzyme (TACE) and a thiol-containing aryl sulfonamide compound (AS). TACE was chosen due to the planar geometry of the active site cleft and its high relevance to the pharmaceutical industry. Due to its role in the release of membrane-anchored cytokines like the tumor necrosis factor-$\alpha$, it is a promising drug target for the treatment of certain types of cancer, Crohn’s disease and rheumatoid arthritis Murphy (2008); Moss et al. (2008b); Rose-John (2013). The ligand under consideration is part of a series of thiol-containing aryl sulfonamides which exhibit potent inhibition of TACE, and is supported by a crystallographic structure Rao et al. (2007). This complex provides an important testbed to benchmark our GBS-enhanced method. As we will show, our method is able to find the correct binding pose without requiring all-atom representation or simulation of the ligand/receptor complex.

The binding interaction graph for the TACE-AS complex is constructed by first extracting all the pharmacophore points on ligand and receptor using the software package rdkit Landrum et al. (2006). To simplify numerical simulations, we identity the relavant pairs of pharmacophore points on the ligand and receptor that are within a distance of $4$Å of each other, and whose label pairs are either hydrogen donor/acceptor, hydrophobe/hydrophobe, negative/positive charge, aromatic/aromatic. After this procedure, we get 4 points on the ligand and 6 points on the receptor and create two labelled distance graphs as illustrated in Fig. 1. The knowledge-based potential is derived by combining information from PDBbind Liu et al. (2014); Wang et al. (2005, 2004), a curated dataset of protein-ligand interactions, and the Drugscore potential Gohlke et al. (2000a, b); Mooij and Verdonk (2005). More details are presented in Appendix C, where the resulting knowledge-based potential is shown in Table S1.

Using this knowledge-based potential, we combine the two labelled distance graphs into the TACE-AS binding interaction graph as shown in Fig. 2. A summary of our graph-based molecular docking approach is shown in Fig. 4, which includes a molecular rendering of the predicted binding interactions of the AS ligand in the TACE binding site using the crystallographic structure of this complex (PDB: 2OI0) LABEL:rao2007novel. These interactions correspond to the maximum vertex-weighted clique in the TACE-AS graph. This set of pharmacophore interactions can be used as constraints in a subsequent round of molecular docking to deduce three-dimensional structures of the ligand-receptor complex LABEL:Hindle2002,verdonk2004. We now study the search of the maximum weighted clique on the TACE-AS graph via a hierarchy of algorithms in increasing order of sophistication. As discussed previously, these are:

Random search: Generate subgraphs at random and pick the cliques with the largest weight among the outputs. 2. 2.

Greedy shrinking: Generate a large random subgraph and remove vertices until a clique is obtained. Vertices are removed by taking into account both their degree and their weight. 3. 3.

Shrinking + local search: Use the output of the greedy shrinking algorithm as the input to a local search algorithm.

These form a hierarchy in the sense that random search is a subroutine of greedy shrinking, which is itself a subroutine of shrinking + local search. For each of these algorithms we compare the performance of standard classical strategies with their quantum-classical hybrid versions introduced in Sec. III.3, where the random subgraph is sampled via GBS. For a fair comparison with GBS-based approaches, the classical data is generated as follows: we first sample a subgraph size $N$ from a normal distribution with the same mean $\langle N\rangle$ and variance $\Delta N^{2}$ as the GBS distribution, then uniformly generate a random subgraph with size $N$.

We begin our analysis with a pure GBS random search. We consider GBS with threshold detectors, which register measurement outcomes as either ‘no-click’ (absence of photons) or ‘click’ (presence of one or more photons). We employ either a brute force approach to calculate the resulting probability distribution or, when that becomes infeasible, the exact sampling algorithm discussed in Refs. Quesada et al. (2018); Gupt et al. (2018). Given the complexity of simulating GBS with classical computers, for simplicity in numerical benchmarking, we first consider the simpler case where the maximum clique size is known, so we can post-select GBS data to have a fixed number of detection clicks. This drastically simplifies numerical simulations (see Appendix A.4 for details), at the expense of disregarding data that would otherwise be present in an experimental setting.

For the TACE-AS binding interaction graph, the largest and heaviest cliques both have eight vertices, so we fix $N=8$. There are a total of 19 cliques of this size in the graph (see also Fig. S1 in the Appendix D). In Fig. 5 we show the outcomes of a numerical experiment where a GBS device has been programmed to sample from the Hafnian of $\Omega A\Omega$, with $\Omega$ as in Eq. (5). For simplicity, we choose $\alpha=1$ in Eq. (5), although performance can be slightly improved with optimized values of $\alpha$. On the other hand, the parameter $c$ does not play any role in the post-selected data, but it does change the overall probability of getting samples of size $N=8$. For comparison, we have also studied a purely classical random search, where each data is a uniform random subgraph with $N$ vertices. We observe only three cliques over $10^{5}$ samples. On the other hand, as shown in Fig. 5, GBS is able to produce roughly $300$ cliques directly from sampling, without any classical post-processing. This indicates that the GBS distribution is indeed favouring cliques with large weights, as intended.

Post-selecting on the number of detector clicks is an unwise strategy when employing real GBS devices because it disregards otherwise useful samples. Moreover, the size of the maximum weighted clique is generally unknown. Instead, we can generate cliques from every sample by employing the shrinking strategy discussed in Section III.3.

In Fig. 6 we study the performance of greedy shrinking with GBS data. These data consist of $10^{4}$ samples obtained from an exact numerical sampling algorithm Gupt et al. (2018). Each sample corresponds to a subgraph and, unlike Fig. 5, here any subgraph size is considered. These results show that with GBS and greedy shrinking – a simple classical post-processing heuristic – it is possible to obtain the maximum weighted clique with sufficiently high probability. Indeed, the histogram in Fig. 6 has a sharp peak corresponding to the clique of maximum size $N=8$ and maximum weight $\approx 3.99$. The success rate in sampling from the max weighted clique is $\approx 12\%$ and the overall sampling rate for $N=8$ cliques is $\approx 19\%$. Greedy Shrinking with purely classical random data is shown in the Supplementary Fig. S2. Although the classical distribution is chosen to have the same mean and variance as the GBS distribution, its performance is considerably worse: the maximum weighted clique is obtained only $1\%$ of the time, compared to $12\%$ for GBS. This shows that GBS with greedy shrinking is already able to find the maximum weight clique of the graph after only a few repetitions.

Finally, we study how the cliques obtained from GBS with greedy shrinking can be enlarged or improved via local search. Fig. 7 shows the performance of the hybrid GBS shrinking + local search algorithm, compared to a classical strategy. The results indicate that GBS not only provides better initial estimates after greedy shrinking (zero iteration steps), but it maintains a significant margin compared to classical strategies as the number of steps is increased. After $k=8$ local expansion steps, the probability of finding the maximum weighted clique is as high as $60\%$, while the classical strategy has a considerably smaller success rate of $<30\%$. After many steps, the success rate saturates: using GBS the success rate gets close to $70\%$, while for the purely classical approach it remains under approximately $35\%$.

The role of noise and squeezing is discussed in Appendix D, where we show that GBS success rate is not diminished by the effect of noise, provided that the amount of squeezing is increased accordingly. Therefore, GBS shrinking and its variant with local search are robust against noise, maintaining a significant margin compared to purely classical strategies.

V Conclusions

We have shown that Gaussian Boson Sampling (GBS) can be employed to predict accurate molecular docking configurations, a central problem in pharmaceutical research and development. This is achieved by first mapping the docking problem to the task of finding large cliques in a vertex-weighted graph, then programming the GBS device to sample these cliques with high probability. This constitutes an example of the viability of near-term quantum photonic devices to tackle problems of practical interest.

Established algorithms for obtaining molecular docking configurations exist, but a challenge arises in the context of industrial drug design where large numbers of candidate molecules must be screened against a drug target. In this case, a fast method for predicting docking configurations is required. In principle, photonic devices such as Gaussian Boson Samplers can operate at very high rates, and may potentially provide solutions in shorter timeframes. Additionally, by sampling better random subgraphs, GBS serves as a technique to enhance the performance of classical algorithms because it increases the success rate of identifying large weighted cliques. This property is relevant and applicable in any context where identifying clusters in graphs is important, beyond applications in molecular docking.

More broadly, our results establish a connection between seemingly disparate physical systems: the statistical properties of photons interacting in a linear-optical network can encode information about the spatial configuration of molecules when they combine to form larger complexes. In other words, we have found that when the interaction between fundamental particles is carefully engineered, they acquire collective properties that can be probed to perform useful tasks. A complete understanding of the capabilities of emerging quantum technologies may thus require further exploration of systems that, even if incapable of universal quantum computation, can still be programmed to exhibit properties that can be harnessed for practical applications.

Acknowledgments

We thank Thomas R. Bromley, Nathan Killoran, Nicolás Quesada, Christian Weedbrook, and Maria Schuld for fruitful discussions.

Appendix A Technical details

A.1 Pure-state Gaussian Boson Sampling

Consider $M$ single-mode squeezed states with squeezing parameter $r_{j}$ injected into an $M\times M$ interferometer described by the unitary matrix $U$. For a pure Gaussian state, the matrix $\mathcal{A}$ entering in Eq. (1) can be decomposed as $\mathcal{A}=B\oplus B^{*}$, where $B=U\bigoplus_{j=1}^{M}\tanh(r_{j})U^{T}$ Hamilton et al. (2017). Therefore, the GBS device can be programmed with any matrix $B$ whose spectrum is contained in $[0,1]$. The adjacency matrix $A$ of a graph normally does not have this spectral property. However, it can always be rescaled as $A=cB+d\operatorname{\mathbbm{1}}$ where $B$ has the desired spectrum and $c,d$ are suitable rescaling constants. Using properties of the Hafnian, namely ${\rm Haf}(A)=c^{N}{\rm Haf}(B)$ and ${\rm Haf}(B\oplus B^{*})=|{\rm Haf}(B)|^{2}$, it was found that a GBS device can be programmed to sample from a distribution $p(S)\propto\frac{|{\rm Haf}(A_{S})|^{2}}{c^{N}}$ Brádler et al. (2018).

A.2 Weighted Gaussian Boson Sampling

We now introduce some technical results to show the properties of weighted Gaussian Boson Sampling.

Lemma 1*.*

Let $B$ be defined as in Eq. (3), where $\Omega$ is diagonal with non-zero diagonal elements $c\omega_{j}$ for some $c$. Then, $B$ has spectrum in $[0,c]$ if $2\max_{j}d_{j}\omega_{j}^{2}\leq c^{-1}$.

Proof.   The proof closely follows that of Lemma 1.7 in Chung and Graham (1997). Since the Laplacian $D-A$ is positive semidefinite, so is $B$. Let $f$ be any vector, then

[TABLE]

where $g=D^{1/2}\Omega f$ and where we used that the spectrum of the normalized Laplacian $\tilde{L}$ is contained in $[0,2]$. ∎

The main result of this section is the following decoupling theorem

Theorem 1*.*

Let $B$ be the matrix defined in (3). Then

[TABLE]

Proof.   The proof is based on the following expansion

[TABLE]

where in (a) we use the definition of the Hafnian where PMP is the set of perfect matchings. Equality (b) follows from the definition of $B$, being the Hafnian of a matrix independent of its diagonal elements. In (b) each $M$ contributes to the sum only if $A_{ij}\neq 0$ for all $(ij)\in M$. In the latter case the product $\prod_{(ij)\in M}(\omega_{i}\omega_{j})$ is a product over all possible $\omega_{j}$, as each vertex is visited in $M$ only one time. In (c) we use the fact that the latter product is independent of $M$. ∎

A.3 Biasing the number of detections

We discuss the role of parameter $c$ in biasing the average output size. Consider a single-mode state with squeezing parameter $r$. Being pure, the $\mathcal{A}$ matrix is written as $\mathcal{A}=B\oplus B^{*}$ and, for a single mode, $B=\tanh(r)$. For maximum squeezing $r_{\rm max}$ we find that $B$ can take any value in $[0,c]$ with $c=\tanh(r_{\rm max})$. The resulting average photon number is then $\langle N\rangle=\sinh(r)^{2}=\frac{c^{2}}{1-c^{2}}$ and the variance is $\Delta N^{2}\propto{\langle N\rangle(1+\langle N\rangle)}$. For multiple modes the expressions are similar, though $B$ is a matrix and $\langle N\rangle=\operatorname{Tr}[\frac{B^{2}}{\mathbb{I}-B^{2}}]$, so the normalization factor can be tuned to provide a higher rate to subgraphs of different sizes $N$. Although the maximum clique size is not known a priori, an estimate, e.g. based on random graphs Bollobás and Erdös (1976), is normally enough as the large variance $\Delta N^{2}\approx{\langle N\rangle(1+\langle N\rangle)}$ assures that different sizes are sampled with sufficiently high rate.

Gaussian boson sampling using click detectors yields a discrete probability distribution over subsets $S_{N}$ of $\{1,\dots,M\}$ of dimension $N$. We write $i\in S_{N}$ if the $i$th detector “clicks” and $i\notin S_{N}$ otherwise. The resulting probability distribution is Quesada et al. (2018)

[TABLE]

where $P_{0}^{i}=|0_{i}\rangle\langle{0_{i}}|$ is the projection into the zero photon state and $P_{1}^{i}=\openone-P_{0}^{i}$. The average number of clicks $N$ is then

[TABLE]

where $\rho_{j}$ is the reduced state on mode $j$. Using the fidelity formula for Gaussian states Banchi et al. (2015) we then get

[TABLE]

where $\sigma_{j}$ is the reduced ($2\times 2$) covariance matrix for mode $j$. The above equation can be solved to bias the number of clicks. When the covariance matrix $\sigma$ depends on the normalization factor $c$, we can use a simple line search algorithms to tune $c$ such that $N[\sigma(c)]$ is equal to the desired value.

A.4 Post-selection

Sampling from Eq. (14) requires the calculation of all $p(S_{N})$ for $N=1,\dots,M$. There are exponentially many of these probabilities $\mathcal{O}(2^{M})$. However, if we are interested in samples of a fixed size $N$, then the number of $p(S_{N})$ with fixed $N$ is $\mathcal{O}(\binom{M}{N})\approx\mathcal{O}(M^{N})$. Each probability requires the evaluation of $\mathcal{O}(2^{N})$ determinants, so the complexity is still exponential as a function of $N$ Quesada et al. (2018). However, focusing on postselection with a certain size $N$ reduces the complexity of brute force approaches from exponential to polynomial, although the degree of this polynomial increases with $N$.

A.5 Selecting parameter $\alpha$

For a complete graph with $2n$ vertices the Hafnian is $h_{n}=\frac{2n!}{n!2^{n}}$. The largest Hafnian for non-complete graphs is obtained by removing an edge from the complete graph. The Hafnian is then $\frac{2n-2}{2n-1}h_{n}$, so this non-optimal graph is penalized by a factor $\frac{2n-2}{2n-1}\simeq 1-\frac{1}{2n}$. A possible choice for $\alpha$ is to avoid a counterbalance of this term, so $1+\alpha w_{\rm tot}<\frac{2n-1}{2n-2}$. Nonetheless, we have numerically observed that, at least for sparse graphs, the parameter $\alpha$ does not have to be carefully chosen, and different values of $\alpha$ provide the expected enhancement for the max weighted clique problem.

Appendix B Graph representations of molecular interactions

In this section, we use $\mathbb{L}$ to denote the set of all labels corresponding to the individual pharmacophore point types and $\kappa$ to denote the potential function $\kappa:\mathbb{L}\times\mathbb{L}\rightarrow\mathbb{R}$ that assigns an interaction strength to each pair of labels from $\mathbb{L}$.

B.1 Labeled distance graph

Definition B.1**.**

Labeled distance graph.

Let $S$ be a set of points in three dimensional space S=\big{\{}(x_{i},y_{i},z_{i})\,|\,i\in I\big{\}} for a given index set $I$ heuristically selecting pharmacophore points of a component (either ligand or the binding site) involved in the binding complex. Then labeled distance graph $G_{S}$ is defined as $G_{S}=(V_{S},E_{S},\omega_{S},\alpha_{S})$ where

[TABLE]

is the set of vertices,

[TABLE]

is the set of edges,

[TABLE]

is the weighting function of the edges and

[TABLE]

is a function assigning a pharmacophore point type to each vertex.

Remark*.*

Any labeled distance-graph is a complete graph with $I$ vertices.

B.2 Binding interaction graph

Let $G_{L}=(V_{L},E_{L},\omega_{L},\alpha_{L})$ be a ligand labeled distance graph and $G_{B}=(V_{B},E_{B},\omega_{B},\alpha_{B})$ labeled distance graph for the binding site. Any pair of vertices $(l_{i},b_{i})\in V_{L}\times V_{B}$ is then called a contact between $G_{B}$ and $G_{L}$.

Definition B.2**.**

$\tau$* flexible contact pair.*

Let $c_{i}=(l_{i},b_{i})$ and $c_{j}=(l_{j},b_{j})$ be contacts between $G_{L}=(V_{L},E_{L},\omega_{L},\alpha_{L})$ and $G_{B}=(V_{B},E_{B},\omega_{B},\alpha_{B})$. Then $(c_{i},c_{j})$ is a $\tau$* flexible* contact pair between $G_{L}$ and $G_{B}$ if and only if $|\,\omega_{L}(l_{i},l_{j})-\omega_{B}(b_{i},b_{j})\,|\leq\tau+2\epsilon$, where $\tau$ is the flexibility constant and $\epsilon$ is the interaction cutoff distance.

Remark*.*

Mutual $\tau$ flexibility of contact pairs is a reflexive and symmetric relation, but not necessarily transitive.

For multiple contact pairs to be realized in the binding pose they have to not violate each other’s geometric constraints and hence be pairwise $\tau$ flexible . In the following graph representation, this corresponds to a clique:

Definition B.3**.**

Binding interaction graph.

Let $G_{L}=(V_{L},E_{L},\omega_{L},\alpha_{L})$ be a labeled distance-graph for a given ligand and $G_{B}=(V_{B},E_{B},\omega_{B},\alpha_{B})$ a labeled distance-graph for a given binding site. The corresponding binding interaction graph $I_{L,B}$ is defined as

[TABLE]

where vertex set $\mathcal{V}$ is the set of the pairs over vertex-sets of $G_{L}$ and $G_{B}$

[TABLE]

and $\tau,\epsilon\in\mathbb{R^{+}}$ are the flexibility threshold constant and interaction cutoff distance. Then

[TABLE]

is a maximal set of $\tau$ flexible contact pairs between $G_{L}$ and $G_{B}$ and $\Omega:\mathcal{V}\rightarrow\mathbb{R}$ is an vertex-weighting function defined as

[TABLE]

which encodes the interaction strength between pharmacophore points corresponding to vertices $v_{l}$ and $v_{b}$.

Remark*.*

Most favourable binding pose of ligand described by labeled distance graph $G_{L}$ and binding site described by labeled distance graph $G_{B}$ corresponds to the heaviest vertex-weighted clique of binding interaction graph $I_{L,B}$.

Appendix C TACE-AS complex

The binding interaction graph used in numerical simulations is constructed as follows. First, using the software package rdkit Landrum et al. (2006) all the pharmacophore points on ligand and receptor are extracted. This results in 11 pharmacophore points on the AS ligand and 243 points on the TACE receptor. These sizes are not too large for future GBS devices, but the classical simulation of GBS is intractable for such large problem instances Gupt et al. (2018). Therefore, to enable numerical simulations, we subselect pharmacophore points based on the true binding pose of AS and TACE according to the following two criteria:

Select pairs of pharmacophore points on ligand and receptor that are within $4$Å distance of each other. 2. 2.

From these pairs select the ones whose label pairs are either hydrogen donor/acceptor, hydrophobe/hydrophobe, negative/positive charge, aromatic/aromatic.

Note that in a realistic scenario the true binding pose would be unknown. However, a similar set of points could be obtained based on knowledge commonly employed in drug discovery. For example, ligand pharmacophore points can be heuristically selected and prior knowledge of the binding site location drastically reduces the number of pharmacophore points on the receptor. To reduce the number of receptor pharmacophore points even further, one could use a sliding window to study different sections of the binding site in isolation. Nevertheless, these reduction techniques will be unnecessary when physical GBS devices are built with enough modes. In the case of the TACE-AS complex, we subselect 4 points on the ligand and 6 points on the receptor and create two labelled distance graphs as illustrated in Fig. 1.

Using PDBbind Liu et al. (2014); Wang et al. (2005, 2004), a curated dataset of protein-ligand interactions, we derive a knowledge-based pharmacophore potential. For all protein-ligand interactions, the pharmacophores on ligand and binding site are extracted with rdkit Landrum et al. (2006) and all pairwise distances are accumulated in a single histogram. Subsequently, the Drugscore potential Gohlke et al. (2000a, b) for the six pharmacophore types (negative/positive charge, hydrogen donor/acceptor, hydrophobe and aromatic ring) is computed from the histogram as outlined in Ref. Mooij and Verdonk (2005). The Drugscore potential yields values in the interval $[0,1]$ whereby favourable interactions are close to [math]. Since we want to encode the correct binding pose in a maximum weighted clique, we reflect the resulting potential,

[TABLE]

such that large values in the potential encode desirable interactions.

The resulting knowledge-based potential is shown in Table S1.

Appendix D Supplementary figures

In Fig. S1 we show the position of all maximum cliques (of size $N=8$) inside the graph, ordered from lightest to heaviest total weight. The figure shows that there are two main clusters in the graph: the top right cluster, generally with light weights, and the bottom left cluster with heavy weights. There are also a couple of intermediate cliques where these two clusters are mixed. The maximum weighted clique is shown in the bottom graph, where from node diameter we observe that it is composed by a heavy six-vertex core and two light vertices. Comparison with Fig. 5 shows that all lightweight cliques have a low occurrence rate in a carefully programmed GBS device.

In Fig. S2 we show the output of Greedy Shrinking with purely classical random data. For a fair comparison with the GBS-based approach shown in Fig. 6, the classical data are generated as follows: we first sample a subgraph size $N$ from a normal distribution with the same mean $\langle N\rangle$ and variance $\Delta N^{2}$ as the GBS distribution, then uniformly generate a random subgraph with size $N$. Although the resulting distribution has the same mean and variance as the GBS distribution, by comparing Fig. S2 and Fig. 6, we see that its performance is considerably worse: the maximum weighted clique is obtained only $1\%$ of the time, compared to $12\%$ for GBS.

In Fig. S3 we study the effect of noise and squeezing. The value $r_{\rm max}=0.9702$ corresponds to an average number of detector clicks $\langle N\rangle\simeq 8$. In the lossy case, for a fair comparison, we have increased the squeezing to $r_{\rm max}=0.9780$ in order to maintain the same average $\langle N\rangle\simeq 8$ and have, accordingly, samples of the same average size. As Fig. S3 shows, the success rate is not diminished by the effect of noise, provided that the amount of squeezing is increased accordingly. As a matter of fact, the noisy version with larger squeezing displays a similar success rate after greedy shrinking (iteration 0). As the iterations increase, the success rate of both noisy and noiseless GBS maintain a significant margin compared to the purely classical strategy. The slightly better performance of the noisy case is due to the larger squeezing that changes the shape of the photon distribution, while keeping comparable photon averages with the noiseless case. This analysis shows that both GBS shrinking and its variant with local search are robust against noise, maintaining a significant margin compared to purely classical strategies.