Boltzmann Sampling by Degenerate Optical Parametric Oscillator Network for Structure-Based Virtual Screening

A structure-based lead optimization procedure is an essential step to finding appropriate ligand molecules binding to a target protein structure in order to identify drug candidates. This procedure takes a known structure of a protein-ligand complex as input, and structurally similar compounds with the query ligand are designed in consideration with all possible combinations of atomic species. This task is, however, computationally hard since such combinatorial optimization problems belong to the non-deterministic nonpolynomial-time hard (NP-hard) class. In this paper, we propose the structure-based lead generation and optimization procedures by a degenerate optical parametric oscillator (DOPO) network. Results of numerical simulation demonstrate that the DOPO network efficiently identifies a set of appropriate ligand molecules according to the Boltzmann sampling law.


Introduction
Structure-based drug design (SBDD) is a powerful tool for identifying lead compounds which can be further optimized or filtered to generate drug candidates [1][2][3].When details of a 3D structure of a target protein are known, SBDD can identify suitable compounds, which have high binding affinity and selectivity with a target protein.Several structure-based virtual screening methods based on the lock-and-key model for protein-ligand interactions have succeeded in identifying hit compounds [4][5][6].A lead optimization procedure from given fixed geometry has also been proposed [7,8].In this method, all possible combinations of atomic species to given geometry are considered and ranked based on the score function which is an interaction energy between protein and ligand molecules.Since this procedure is a notorious combinatorial optimization problem, an exhaustive search is used to rank the combinations of atomic species in order to discover drug candidates.
Recently, there has been several attempts to solve computationally hard combinatorial optimization problems by using real physical systems rather than a conventional digital computer [9][10][11].Such systems take the approach of physical simulation of Ising models, to which various combinatorial optimization problems can be mapped [12][13][14].These systems are expected to solve combinatorial optimization problems as a random sampling machine like Markov chain Mote Carlo (MCMC) simulation which is often used for machine learning [15][16][17].The degenerate optical parametric oscillator (DOPO) network is one of such systems, which is recently proposed and demonstrated as a coherent Ising machine (CIM) [18][19][20].The advantage of CIM is that any graph topology, including all to all connections, can be implemented without a time-consuming task of graph transformation.The performance of DOPO-based CIM is evaluated against NP-hard MAX-CUT problems which are equivalent to find the ground state of the Ising Hamiltonian without Zeeman term [18,21].However, applications of CIM to real-world problems have not been proposed yet.Here, we show the mapping of the structure-based lead optimization problem onto CIM and present the benchmark simulation results by the DOPO network.

Computation Principle of CIM
CIM has been proposed to solve combinatorial optimization problems by exploiting the phase transition of a network of quantum oscillators.This machine searches the spin configuration {σ i } which corresponds to the ground state of the Ising Hamiltonian H representing cost function of an optimization problem.The Isnig Hamiltonian H is defined as where σ i ∈ {±1} is a binary spin state of site i and interactions between spins are described by the coupling matrix {J ij }.A time-division-multiplexed (TDM) pulsed DOPO with an optical fiber ring cavity can be a leading candidate for a scalable CIM [20].Figure 1 shows the basic configuration of the CIM based on TDM pulsed DOPO [21].The Ising spin σ i is represented by the in-phase amplitude c i of the i-th DOPO pulse.The positive c i value corresponds to up-spin (σ i = +1) while the negative c i value corresponds to down-spin (σ = −1).Each DOPO pulse is in a squeezed vacuum state at a pump rate below oscillation threshold, so that the positive c i (up-spin) and negative c i (down-spin) co-exist as a superposition state.This allows the quantum parallel search for the ground state of the Ising Hamiltonian.The Ising coupling J ij is implemented by optical homodyne detection and electric feedback control.The inferred value ci by the optical homodyne measurement is used to compute an appropriate coupling field to other spin σ j together with the initially memorized J ij values in the digital circuit.An actual feedback pulse for the j-th pulse is produced by the optical amplitude/phase modulation based on the feedback circuit output for a part of the pump laser pulse which is phase-coherent with the DOPO pulse due to the phase sensitive amplification of DOPO.
The dynamics of the in-phase component c and quadrature-phase component s of the field of the i-th DOPO pulse is described as [22,23] where ξ ij is the coupling constant from the j-th to the i-th DOPO pulse and p is the pump rate.
where σ i = sgn(c i ), N is the number of DOPO pulses and is the order of ξ ij .If we ignore the higher order correction given by the third term of Equation ( 4) and exclude the constant term N, the overall photon decay rate Γ is identical to the Ising Hamiltonian.Therefore, an oscillation mode {σ i } which achieves the minimum overall photon decay rate Γ provides the ground state of the Ising Hamiltonian.

Implementation of Zeeman Terms
The Ising Hamiltonian with a Zeeman term induced by an external magnetic field is represented as where h i is the local external field to spin i.In a wide variety of NP problems, they can be mapped to the Ising model with Zeeman terms [12].Therefore we need to expand the argument of the previous subsection to the Ising model with Zeeman terms.This can be done by adding an extra term to the Equation (2) as where |c| = (1/N) ∑ i |c i | is the average absolute amplitude overall DOPO pulses at time t and h is a constant that determines the strength of the Zeeman term.Here, we ignore the dynamics of the quadrature-phase component due to the reason mentioned above.From this equation, the overall photon decay rate is computed, which is equivalent to Equation ( 5) when h equals 1.

Lead Optimization Procedure
The lead optimization procedure used in this work is based on the paper by Ogata et al. [8].It is assumed that the coordinates of the compound structures, which are called geometry here, are extracted from a high resolution structure imaging of the protein-ligand complex.Then all the coordinates in the geometry are classified according to their bond order types (sp3, sp2, etc.) and atomic species (CH3, CH2, CH, NH2, NH, etc.).For instance, the geometry determined by a target protein is given as follows: where X, Y, and Z represent the coordinates in the geometry and "• • • " is a generic representation of the bonds which connects the atomic species.Replacing Y with an atomic species " =CH−" may generate a chemically incomplete compound X=CH−Z for which "=" and "−" indicate a double and a single bonds, respectively.Then, X must be assigned to an atomic species connected through a double bond to Y (e.g., O= or CH2=).Similarly, Z must be assigned to an atomic species capable of linking to Y through a single bond (e.g., −CH3 or −NH2).By assembling all possible combinations of these atomic species, four compounds are obtained: O=CH−CH3, O=CH−NH2, CH2=CH−CH3 and CH2=CH−NH2.
After assigning atomic species to the geometry, the interaction energy between the obtained compounds and the target protein must be evaluated.The cost function, including an interaction energy and solvation energy, is calculated for each atomic species separately as where i is the position to which an atomic species is assigned and n is the atomic species.E v in is the van der Waals interaction energy between the atomic species and the target protein.E e in is the electrostatic interaction energy between the two, E h in is the hydrogenic bond energy and E s in is the solvation energy.The total cost function is represented as E total = ∑ i E in i , where n i is the atomic species which is assigned to the position i.The details of the interaction energy are described in the report by Ogata et al. [8].In the original method, this cost function is calculated for all possible combinations of atomic species and all such ligand compounds are ranked in the ascending order according to the cost function value.Since this approach experiences combinatorial explosion as the problem size becomes large, we take the approach to sample only ligand compounds with low cost function values from all possible combinations.Note that there are ∼3 × 10 12 ligand compounds even for such a small problem size with N position = 16 and N atom = 6 atomic species where N position and N atom are the numbers of sites and atomic species, respectively.

Mapping to the Ising Hamiltonian
The problem considered above is summarized in the following constrained optimization problems: From given geometry and atomic species, sample the low energy assignments {i, n} of atomic species to geometry satisfying the constraints of the bond order.We can map the problem to the Ising Hamiltonian with Zeeman terms as follows.
Let us consider a binary variable x in ∈ {0, 1} with position i on the geometry and atomic species n: x in = 1 if atomic species n is assigned to position i, When the total number of atomic species and positions are N atom and N position , respectively, the number of variables x in required to compute the cost function is N = N atom N position .The penalty between x in and x jm is now defined as J in;jm = 1 if position i is adjacent to j and the bond order between n and m are inconsistent, 0 otherwise.
The total cost function with x in and J in;jm for the considered problem is described by where A, B and C are adjustable constants which determine the contribution weights of each term in the cost function and E in is defined by Equation ( 7).The fist term in the equation comes from the constraint that exactly one atomic species must be assigned to each position.The second term denotes the constraint of the bond order.The first and second terms must be zero for a legitimate compound satisfying the two constraints.Finally, the third term is the interaction energy which we want to minimize under the two constraints.We appropriately choose the parameters A, B and C, to efficiently find the ground states on low-energy states of the Hamiltonian under the condition of satisfying the constraints.Since the above formulation (Equation ( 10)) uses the binary variable x in ∈ {0, 1}, we need to transform the variable x in to σ in ∈ {−1, 1} by the relation x in = (σ in + 1)/2 for mapping the cost function on the CIM.Equation ( 10) can be then represented by where Since the constant term doesn't change the energy landscape, we can neglect the third term of Equation (11).The resulting cost function is the Ising Hamiltonian with Zeeman terms, so that we can implement it directly on the CIM.

Several Heuristic Modifications
A CIM based on the DOPO network simulates an Ising model and reports not only the ground states but also low-energy excited states.DOPO pulses, which represent Ising spins σ i ∈ {−1, 1}, can actually take continuous values as the amplitude of the in-phase component c, whereas the original Ising spins take two discrete spin states.It has been found that CIM solves MAX-CUT problems without a Zeeman term with high performance [18,21].However, our problem has a Zeeman term and the goal is not only finding the ground states with an exact minimum energy, but rather sampling many low-energy states which can stably bind to the target protein and thus become drug candidates.
Here, we introduce three heuristic methods to improve the performance of the CIM to solve the given problems.A first heuristic method is the gradual pumping where the pump parameter p is gradually increased from below to above the oscillation threshold.The pump parameter p acts as a gain, and when there are no mutual couplings, the amplitude of each DOPO pulse at the steady state is proportional to p − 1 for p > 1 .The gradual pumping makes quantum noise effects gradually weaker as the amplitude grows .In this paper, pump p is increased linearly with time until some limit value.
As a second heuristic method, we replace the amplitude of the coupling field in Equation ( 6) as By this replacement, each DOPO pulse experiences the effect of interaction as if all other DOPO pulses are oscillating with the same amplitude.This modification allows us to avoid the breakdown of the mapping protocol due to the inhomogeneous DOPO amplitudes.Finally, we consider the relaxation of {J in;jm } (Equation ( 9)) to a continuous value.When the CIM solves an Ising problem with a Zeeman term, we add it with a constant h as described in Equation ( 6).The parameter h empirically tunes the balance between the Ising interaction term and the Zeeman term.If there are large variation in Zeeman coefficients h in , the Zeeman Hamiltonian dominates over the Ising Hamiltonian and CIM picks up biased solutions on the same Zeeman energy surface.In our problem, the large variation in h in is caused by the term ∑ j,m J in;jm .To avoid this situation, we introduce the relaxation of {J in;jm } to a continuous value while keeping symmetry, zero elements and non-zero elements with positive real values so that the second term of Equation ( 13) satisfies Since J in;jm is only related to the second constraint which should be zero for legitimate states, this relaxation doesn't affect the energy of solutions.If we represent J as the vector whose elements are non-zero elements of {J in;jm }, the problem is finding J which satisfies where 1 is a constant vector whose elements are one, W is an N nonzero × N matrix which specifies the locations of non-zero elements of {J in;jm } and N nonzero is the number of non-zero elements of {J in;jm }.
There are infinite ways of this relaxation because N nonzero is larger than N in general.In this paper, we use J which is found by calculating the Moore-Penrose pseudoinverse of W [24].

Numerical Simulation Results
We tested our proposed protocol including heuristic modifications by the numerical simulation for the simple benchmark problem.The benchmark we used here is a 6-membered ring placing near the artificially generated Ala-Asp-Ala tripeptide (Figure 2a).To generate compounds, six types of atomic species, =CH-, -CH=, -NH-=N-, -N= and -O-, are assigned to the ring.Therefore, the number of spins N which CIM needs to solve the problem is N = N atom N position = 36, where N atom = N position = 6.The interaction energies for all possible compounds are calculated in advance.Benzene and pyridine are detected as the configuration with the lowest and 2nd lowest energies, respectively, in 0.53 kcal/mol score difference (Figure 2b,c).The details of the interaction energy E in and J in;jm of this benchmark are available in the Supplementary Materials.At first, we tested our constraint Hamiltonian without the score term, i.e., the first and second terms of Equation ( 10) without the interaction energy (C = 0), to investigate the parameter dependence of the probability of satisfying the constraints.Figure 3 shows the success probability of satisfying the constraints over such parameters as pump limit p, the weight of Zeeman term h and the ratio of the bond constraint over the position constraint B/A.The success probability depends strongly on h rather than B/A and p.The optimal value for the weight of the Zeeman term is h ∼ 0.6.Next, we investigate the statistics for the interaction energy, the third term of Equation ( 10), with p = 1, h = 0.6 and B/A = 0.6.We plotted the histograms of the interaction energy for various C/A values only for solutions which satisfy the constraints, given by the first and second terms of Equation (10), in Figure 4.The corresponding Boltzmann distribution at properly chosen temperature is also plotted in red.If we represent the Ising Hamiltonian of Equation ( 10) as H(x) = H cons (x) + CH score (x), the Boltzmann distribution which satisfy the constraints is represented by where E is the interaction energy, β n = C/k B T is the normalized temperature parameter, k B is the Boltzmann constant, T is the temperature and D(E) is the number of states at the energy E. Changing the parameter C corresponds to changing the normalized temperature parameter for the Boltzmann distribution.We estimated the normalized temperature parameter for each C/A by the least squares method.The observed distributions are well matched with the Boltzmann distributions although the success probability of satisfying the constraints decreases as C becomes large.Finally, we plotted the histograms for finding degenerate states at the energy interval E = [−3.854,−3.768) and E = [−4.626,−4.540) in Figure 5 in order to investigate the random sampling capability for the solutions on the same energy surface.All degenerate states are found with almost equal probability at C/A = 0.5.Since our purpose is to sample all possible compound configurations with low interaction energies, we confirmed that the optimal parameter of C/A is ∼0.5 in terms of both low energy temperature and random sampling capabilities.

Conclusions
In this paper, we proposed the lead optimization protocol for the CIM and tested the method for a simple benchmark problem.This is the first benchmark study of the DOPO-based CIM that the problem Hamiltonian includes Zeeman terms.We confirm that the CIM satisfies the constraints with high success probability and the resulting histograms of obtaining various compound configurations are well matched with the Boltzmann distribution with low energy temperature.The normalized temperature parameter of the Boltzmann distribution can be tuned by changing the weight C of the interaction Hamiltonian.
We also investigated randomness of the solutions on the same energy surface.All the solutions in the same energy interval are covered with almost identical probability.Since exact randomness and the Boltzmann distribution are not needed for our purpose, the present results were acceptable.
We do not need the same solutions which are already obtained, so that CIM should pick up different solutions in subsequent runs.There are a lot of similar, further tricks that we can employ to improve the performance of our method [25].

Supplementary Materials:
The data which we used for the benchmark study is available online at www.mdpi.com/1099-4300/18/10/365/s1, Table S1: Interaction Energy E in of benchmark problem, Table S2: J in;jm of benchmark problem.

Figure 1 .
Figure 1.Coherent Ising machine (CIM) based on time-division-multiplexed (TDM) pulsed degenerate optical parametric oscillator (DOPO) with measurement and feedback control.Both local oscillator (LO) pulses and feedback (FB) pulses are taken from the pump laser.A parametric gain is provided by a periodically-poled LiNbO 3 (PPLN) waveguide device and an optical ring cavity is formed by a fiber with ∼1 km length.

Figure 3 .
Figure 3.The success probability of satisfying the constraints for 1000 identical trials.The parameter p is the final pump rate for gradual pumping.The parameters of the Hamiltonian (Equation 10) are set to A = 1 and C = 0 for all the results.

Figure 4 .
Figure 4. Histograms of the interaction energies of the final states of CIM over 1000 runs.The blue bar is the simulation result and the red dot is the estimated Boltzmann distribution.The at the bottom of each figure shows the number of for the given Hamiltonian.All histograms are normalized to 1.
The last terms on the right hand side of the equations describe pump quantum noise and vacuum fluctuation incident from the open port of the output coupler.A s is the amplitude of DOPO at a normalized pump rate p = 2, where p = 1 corresponds to the threshold pump rate.dW c i and dW s i are independent zero-mean Gaussian noise processes.It is clear from Equation (3) that the quadrature-phase component of an each DOPO pulse decays exponentially and only the in-phase component builds up.At the steady state, the total saturated gain of the network ∑ i (p − c 2 i − s 2 i ) equals the overall photon decay rate Γ.If we ignore noise terms and quadrature-phase components, and assume |ξ ij | |c [18] it can be shown analytically that Γ is described as[18]