Paper deep dive
Molecular Design beyond Training Data with Novel Extended Objective Functionals of Generative AI Models Driven by Quantum Annealing Computer
Hayato Kunugi, Mohsen Rahmani, Yosuke Iyama, Yutaro Hirono, Akira Suma, Matthew Woolway, Vladimir Vargas-Calderón, William Kim, Kevin Chern, Mohammad Amin, Masaru Tateno
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 95%
Last extracted: 7/21/2026, 3:37:21 AM
Summary
This paper introduces a novel framework for molecular design using deep generative models integrated with D-Wave quantum annealing. The authors propose a Neural Hash Function (NHF) that serves as both a regularization and binarization scheme, enabling the transformation between continuous and discrete signals in the objective function. The resulting quantum-annealing generative models produce small molecules with higher validity and drug-likeness scores than fully classical models and even the original training data, demonstrating the potential of quantum annealing to enhance feature space sampling in drug discovery.
Entities (8)
Relation Signals (6)
D-Wave → implements → Quantum Annealing
confidence 98% · integrating with a D-Wave annealing quantum computer... Quantum annealing... implement these dynamics in hardware
Neural Hash Function → usedin → Quantum Annealing
confidence 95% · NHF is used both as the regularization and binarization schemes... in the error evaluation (i.e., objective) function
Quantum Boltzmann Machine → usedaspriorfor → Variational Autoencoder
confidence 94% · utilized QBM as a prior distribution of the compound generation model
Neural Hash Function → enables → binarization
confidence 93% · NHF... used both as the regularization and binarization schemes simultaneously
Quantum Annealing → improves → drug-likeness
confidence 92% · compounds generated via the quantum-annealing generative models exhibited higher quality in both validity and drug-likeness
ChEMBL → usedtotrain → Variational Autoencoder
confidence 90% · We trained the model using the ChEMBL public dataset
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Deep generative modeling to stochastically design small molecules is an emerging technology for accelerating drug discovery and development. However, one major issue in molecular generative models is their lower frequency of drug-like compounds. To resolve this problem, we developed a novel framework for optimization of deep generative models integrated with a D-Wave quantum annealing computer, where our Neural Hash Function (NHF) presented herein is used both as the regularization and binarization schemes simultaneously, of which the latter is for transformation between continuous and discrete signals of the classical and quantum neural networks, respectively, in the error evaluation (i.e., objective) function. The compounds generated via the quantum-annealing generative models exhibited higher quality in both validity and drug-likeness than those generated via the fully-classical models, and was further indicated to exceed even the training data in terms of drug-likeness features, without any restraints and conditions to deliberately induce such an optimization. These results indicated an advantage of quantum annealing to aim at a stochastic generator integrated with our novel neural network architectures, for the extended performance of feature space sampling and extraction of characteristic features in drug design.
Tags
Links
- Source: https://arxiv.org/abs/2602.15451v3
- Canonical: https://arxiv.org/abs/2602.15451v3
Trouble viewing inline? Open PDF directly →
Full Text
72,491 characters extracted from source content.
Expand or collapse full text
[1] 1] to Implementation Laboratories, Central Pharmaceutical Research Institute, Tobacco Inc., 1-1 Murasaki-cho, , 569-1125, 2] -Wave Systems Inc., 3033 Beta Ave, , 5G 4M9, Molecular Design beyond Training Data with Novel Extended Objective Functionals of Generative AI Models Driven by Quantum Annealing Computer -Calderón tateno1611@gmail.com [ [ Abstract Deep generative modeling to stochastically design small molecules is an emerging technology for accelerating drug discovery and development. However, one major issue in molecular generative models is their lower frequency of drug-like compounds. To solve this problem, we develop a novel, quantum annealing, generative-model approach for optimizing deep generative models, an approach that includes integrating with a D-Wave annealing quantum computer. Of particular note, our neural hash function (NHF) is used simultaneously as a regularization scheme and binarization scheme, of which the latter is for transformation between continuous and discrete signals of the classical and quantum neural networks, respectively, in the error evaluation (i.e., objective) function. The compounds generated via the quantum-annealing generative models exhibit higher quality in both validity and drug-likeness than those generated via the fully-classical models, and even exceed the training data in terms of drug-likeness features, without any restraints and conditions to deliberately induce such an optimization. These results suggest that quantum annealing could be used as a stochastic generator integrated with our novel neural network architectures for extending the performance of feature space sampling and the extraction of characteristic features in drug design. 1 Introduction In the field of drug discovery, efficiently designing molecular structures with optimal chemical properties and synthetic accessibility is a complex and important research area. Typical approaches such as iterating through the design and experiment cycle can only search a small region in a vast chemical space, where the number of synthesizable molecules of size acceptable for drugs is estimated to be over 1060 [bohacek_art_1996]. Machine learning and deep learning-based drug design can efficiently explore a wider region in such a huge chemical space. Actually, by applying recent achievements in deep learning and deep generative models, various deep generative models have been reported for molecules with desired chemical properties [ajagekar_molecular_2023, dollar_attention-based_2021, zhao_comprehensive_2025]. Despite these efforts, two challenges with compounds generated from existing generative models remain: 1) the low frequency of “drug-like” molecules that satisfy the activity to the target proteins and have acceptable chemical properties and synthetic accessibility; and 2) a trade-off between the character of the compounds and the diversity in the chemical structure [lambrinidis_challenges_2018]. One cause of these problems is the lack of data, resulting in low generalization due to overfitting. Currently, the number of compounds that can be synthesized is ∼ 1010 [sadybekov_synthon-based_2022], whereas the expected size of the chemical space is 1060 (as previously mentioned)—thus, the number of synthesizable compounds represents only a very small proportion of the entire chemical search space, and is insufficient for use as training data. Quantum machine learning (QML) (also referred to as quantum artificial intelligence (QAI)) is a growing area of research that combines quantum computing and machine learning. QML investigates how quantum resources such as superposition, entanglement, and tunneling can accelerate or enhance classical learning models. Most early work has focused on the gate-based paradigm, where data and model parameters are encoded into quantum circuits composed of unitary transformations and projective measurements. Algorithms such as quantum support vector machines, kernel estimators, and quantum neural networks are implemented through parameterized quantum circuits, optimized using hybrid quantum-classical feedback loops [devadas_quantum_2025]. However, in many common scenarios, training parameterized quantum circuits has been shown to be unlikely to be able to scale due to trainability issues known as barren plateaus [mcclean_barren_2018]; it may also be possible for classical computing to efficiently simulate such quantum circuits [cerezo_does_2025]. An alternative and more physically motivated framework arises in quantum annealing [albash_adiabatic_2018, johnson_quantum_2011], which provides an analog realization of optimization and sampling tasks central to machine learning. In quantum annealing, learning problems are mapped onto an Ising Hamiltonian and the system searches for the low-energy configurations. Sampling from these configurations can be guided by training the parameters of the Ising Hamiltonian to learn a “binary” data distribution or a “binary” latent representation of an autoencoder. Quantum annealers, such as D-Wave’s Advantage2TM quantum computer, implement these dynamics in hardware, enabling large-scale exploration of complex energy landscapes for machine learning applications [amin_quantum_2018, benedetti_quantum-assisted_2017, winci_path_2020]. Recent experiments have demonstrated that quantum annealing achieves a scaling advantage in approaching low-energy configurations [king_quantum_2023, king_beyond-classical_2025], while the resulting sampling distributions cannot be efficiently reproduced by any known classical simulation [king_beyond-classical_2025]. Therefore, it is important to explore ways to harness this intrinsic sampling power for quantum-enhanced generative modeling in areas such as drug discovery and materials design. In this work, our approach started from Variational Autoencoder (VAE)-based generative models for generation and inference of chemical compounds. VAE sets the approximated posterior distribution of the latent variables and optimizes the evidence lower bound (ELBO) instead of the true log-likelihood, which is generally intractable. By the amortized inference and the reparameterization framework [kingma_auto-encoding_2013, rezende_stochastic_2014], VAE can efficiently train its objective function, which is also referred to as the loss function with the reconstruction loss and regularization terms included, and is widely used for the generative model. TransVAE [dollar_attention-based_2021] also involves the VAE that generates character sequences of molecules via a combination of Transformer-based Encoder/Decoder and continuous-valued latent space. More recently, VAE with discrete latent variables, Discrete VAE (DVAE) [rolfe_discrete_2016], in which the generative process is driven by a Boltzmann Machine (BM), was adopted for the generation of chemical compounds [gircha_hybrid_2023]. DVAE incorporates discrete latent variables with a discrete prior, making it a more suitable alternative for modeling data with categorical structures such as molecular descriptors tokenized by employing a particular transformation scheme. In addition, DVAE is suitable for combining a VAE architecture with quantum computing in which measurement outcomes are binary (spin) states. Quantum VAE [khoshaman_quantum_2018] replaced the prior distribution from a classical BM with a quantum Boltzmann machine (QBM) [amin_quantum_2018] in which the D-Wave quantum annealer leveraged the sampling from the Boltzmann distribution. In DVAE, converting a continuous representation of the input data to a binary representation in the final encoding step, is not differentiable and, thus, prevents gradients of the error (reconstruction loss) from being backpropagated to the encoder. To solve this issue, DVAE introduced auxiliary continuous variables for training the approximating posterior, which stochastically converted from the discrete variables [rolfe_discrete_2016, khoshaman_gumbolt_2018]. Considering a specific form of the smoothing distribution, the error can be propagated back through the latent variables by reparameterization. During training, only discrete variables are used in the prior distribution, whereas continuous variables are used in the approximating posterior. Thus, non-differentiability for backpropagation is still a crucial issue in this scheme. We propose another approach to solve the non-differentiability of the encoder by involving our novel scheme in the objective function of the autoencoder-based generative model. In this scheme, inspired by Deep Hashing [erin_liong_deep_2015], the output of the encoder is binarized using outputs of a “deterministic” function (i.e., a stochastic distribution is not involved) as the latent variables. An additional term that aims to reduce the loss in binarization is included in the total objective (loss) function for error evaluation of the system. Carefully defining the structure of the loss function and its derivatives enables the loss to backpropagate through the latent variables without smoothing and stochastic reparameterization (differentiability). We applied this scheme to our generative models of chemical compounds, thereby showing that employment of this scheme improved the validity and drug-likeness of the generated compounds in both classical and quantum computations of the prior distributions. In addition, we utilized QBM as a prior distribution of the compound generation model, thus indicating that the quantum prior outperforms the classical counterpart in the quality of the generated compounds. Interestingly, the compounds generated via our “quantum” generative model exhibited a higher quality in both structural validity and drug-likeness scores than those of even the training dataset without any conditions to deliberately induce optimal trends of the evaluation scores. Thus, our present analysis suggests that our novel objective functionals integrated with the D-Wave quantum annealer should possess a powerful potential for the sampling task of the drug discovery field. 2 Results 2.1 Quantum Annealing To solve the learning problems using quantum annealing, the Ising Hamiltonian is formulated as H(s)=A(s)HD+B(s)HPHD=−∑iσix,HP=∑i,jJijσizσjz+∑ihiσiz splitH(s)&=A(s)H_D+B(s)H_P\\ H_D&=- _i _i^x, 21.68121ptH_P= _i,jJ_ij _i^z _j^z+ _ih_i _i^z\\ split (1) where HDH_D is a transverse-fiels driver promoting quantum tunneling, HPH_P encodes the cost function or model parameters and σix,σiz _i^x, _i^z are Pauli operators for i-th element. By slowly varying the annealing schedule A(s)A(s) and B(s)B(s), the system evolves toward the low-energy configurations of HPH_P. Sampling from these configurations can be guided by training the parameters JijJ_ij and hih_i to learn the prior distribution of binary latent variables (see Equation 2). 2.2 VAE and DVAE The training framework of the original VAE is formulated by the variational inference, which maximizes ELBO on the true log-likelihood of data: ℒ=qϕ(|)[logpθ(|)]−DKL(qϕ(|)∥pψ())≦logpθ(),L=E_q_φ(z|x)[ p_θ(x|z)]-D_KL (q_φ(z|x)\ \|\ p_ψ(z) ) p_θ(x), (2) where ϕ,θφ,θ and ψ denote the trainable parameters for the encoder, decoder and prior. The prior parameters ψ includes the parameters of Ising Hamiltonian Jij,hiJ_ij,h_i (see Equation 1). The first term corresponds to the reconstruction loss which measures the expected log-likelihood (decoder) pθ(|)p_θ(x|z) of the data x given the latent variables z under the approximate posterior (encoder) qϕ(|)q_φ(z|x). The second term represents the Kullback-Liebler (KL)-divergence between the approximated posterior qϕ(|)q_φ(z|x) and the prior pψ()p_ψ(z). The training objective is argminϕ,θ,ψ∼[Lelbo(;ϕ,θ,ψ)] φ,θ,ψargmin\ \ E_x [L_elbo(x;φ,θ,ψ)] (3) Lelbo≔−qϕ(|)[logpθ(|)]⏟Lrec(;ϕ,θ)+DKL(qϕ(|)∥pψ()) L_elbo -E_q_φ(z|x)[ p_θ(x|z)]_L_rec(x;φ,θ)+D_KL (q_φ(z|x)\ \|\ p_ψ(z) ) (4) where D denotes the data (empirical) distribution. DVAE [rolfe_discrete_2016] incorporates discrete latent variables with a discrete prior, making DVAE a better alternative for modeling data with categorical structures. Simplified Molecular Input Line Entry System (SMILES) [weininger_smiles_1988] and SELF-referencing Embedded Strings (SELFIES) [krenn_selfies_2022] are widely used for encoding molecules to sequences of strings. We created molecular descriptors tokenized from SMILES strings. In addition, DVAE is suitable for combining the VAE architecture with quantum computing where the observed values are binary (spin) states. DVAE proposed a stochastic binarization to characterize the posterior distribution qϕ(|)q_φ(z|x) for D-dimensional discrete latent variables =(zi)i=1,⋯,D∈0,1Dz=(z_i)_i=1,·s,D∈\0,1\^D from outputs of the encoder =(li)i=1,⋯,Dl=(l_i)_i=1,·s,D as follows: qi=σ(li)zi=ℋ(qi+ρi−1)=ℋ(li+σ−1(ρi)),ρi∼Unif(0,1), splitq_i&=σ(l_i)\\ z_i&=H(q_i+ _i-1)=H (l_i+σ^-1( _i) ), 14.45377pt _i (0,1),\\ split (5) where σ(x)=11+e−x,σ−1(y)=logy−log(1−y)σ(x)= 11+e^-x,σ^-1(y)= y- (1-y) denote the sigmoid function as its inverse, and ℋH denotes the Heaviside step function. The main obstacle of VAE with discrete latent variables is that the function z=zϕ(x,ρ)z=z_φ(x,ρ) is non-differentiable such as Heaviside, which prevents backpropagation using the reparameterization trick. Gumbolt [winci_path_2020, khoshaman_gumbolt_2018] applied the Gumbel trick [maddison_concrete_2016]: a smoothing of z to obtain a continuous variable =(ζi)i=1,⋯,D∈(0,1)D ζ=( _i)_i=1,·s,D∈(0,1)^D by ζi=σ(li+σ−1(ρi)τ) _i=σ ( l_i+σ^-1( _i)τ ), where τ is a smoothing parameter with → ζ in the limit of τ→0τ→ 0. 2.3 Involvement of the Neural Hash Function (NHF) Herein, we aimed to solve the afore-mentioned, non-differentiable features derived from the binarization of discrete latent variables. We first constructed a Transformer [vaswani_attention_2017]-based encoder-decoder architecture with discrete latent variables. Let ∈ℕN×CX ^N× C represent the tokenized SMILES of N molecules using a vocaburary size of V, where C denotes the maximum number of tokens, Each element (Xnc∈1,2,⋯,VX_nc∈\1,2,·s,V\) of X represents the index of the c-th token of the n-th molecule. Tokens are embedded in a dmodeld_model-dimensional space and fed into an encoder fϕf_φ including Transformer layers, a Neural Tensor Network block, and MLP layers (see Methods) to obtain the D-dimensional fixed length vectors =(1,2,⋯,N)∈ℝN×DH=(h_1,h_2,·s,h_N) ^N× D, Binary latent variables =(1,2,⋯,N)∈0,1N×DZ=(z_1,z_2,·s,z_N)∈\0,1\^N× D are obtained from H through a binarization function and fed into a decoder gθg_θ including another Neural Tensor Network and Transformer layers. Finally, outputs of the decoder are passed through a softmax layer to get probabilities of reconstructed tokens ^∈ℝN×C×V X ^N× C× V. Inspired by Deep Hashing [erin_liong_deep_2015], we derived the loss function LnhfL_nhf from the ELBO of the joint distribution of the binary latent variables z and the continuous ones h (for details, see the Discussion section), referred to here as the neural hash function (NHF), as the objective (loss) function of our generative models (see Methods for details): Lnhf L_nhf =Lrec(;θ,ϕ)+Lprior(;ϕ,ψ)+Lquant(;ϕ) =L_rec(X;θ,φ)+L_prior(Z;φ,ψ)+L_quant(Z;φ) (6) Lrec(;θ,ϕ) L_rec(X;θ,φ) ≔1N∑n=1NCE(n,^n) 1N _n=1^NCE(X_n, X_n) (7) Lprior(;ϕ,ψ) L_prior(Z;φ,ψ) ≔−1Nlogpψ(n) - 1N p_ψ(z_n) (8) Lquant(;ϕ) L_quant(Z;φ) ≔λfro2N‖−‖F2+λortho2∑l=1L‖llT−‖F2, _fro2N\|Z-H\|_F^2+ _ortho2 _l=1^L\|W_lW_l^T-I\|_F^2, (9) where CECE denotes the cross entropy and ∥⋅∥F\|·\|_F denotes the Frobenius norm (L2L_2 norm) defined for matrices. LrecL_rec is the reconstruction loss as the cross entropy between inputs and reconstructed tokens. LpriorL_prior is the approximated cross entropy of the prior distribution pψ()p_ψ(z) and the empirical distribution. LquantL_quant means the quantization loss to produce good binary codes. The second term encourages the weight matrices l,l=1,2,⋯,LW_l,l=1,2,·s,L of the L-layer MLP in the encoder to be orthogonal, that is, requiring the different dimensions to be independent. In Deep Hashing, the transformation from H to Z is performed by a non-smooth function such as Heaviside function. To optimize the model parameters by the stochastic gradient descent method, we defined the gradient of LnhfL_nhf with respect to the parameters of the decoder θ, prior ψ and encoder ϕφ as follows (see Methods for details): ∂Lnhf∂θ ∂ L_nhf∂θ =∂Lrec(;θ,ϕ)∂θ = ∂ L_rec(X;θ,φ)∂θ (10) ∂Lnhf∂ψ ∂ L_nhf∂ψ =∂Lprior(;ϕ,ψ)∂ψ = ∂ L_prior(Z;φ,ψ)∂ψ (11) ∂Lnhf∂ϕ ∂ L_nhf∂φ =∑n=1NnT∂n∂ϕ = _n=1^N δ_n^T _n∂φ (12) n δ_n ≔λfroN(n−n)+∂Lrec(;θ,ϕ)∂n+∂Lprior(;ϕ,ψ)∂n _froN(z_n-h_n)+ ∂ L_rec(X;θ,φ) _n+ ∂ L_prior(Z;φ,ψ) _n (13) The scheme of the entire model and the workflow of the training and generation phase are shown in Figure 1. Figure 1: Scheme of the Neural Hash Function (NHF)-based generative model. (a) Input compound X is represented by employing SMILES sequences. Each token from SMILES is embedded in a 160-dimensional vector and added the positional encoding. The embedded tensor is fed into the Transformer encoder layers and flattened by the Neural Tensor Network layer (NTN1) to obtain the continuous vector. h is transformed by our NHF to the binary vector, and then transformed into the latent matrix through another NTN layer (NTN2). In decoder block, the input tensors are passed through the self-attention layer with the subsequent masks. The cross-attention between the output of the self-attention layer and the latent matrix. Finally, SMILES sequences are reconstructed from outputs of the decoder through the softmax layer. The reconstruction loss between the input and reconstructed tokens (LrecL_rec), the cross entropy of the prior distribution (LpriorL_prior) and the quantization loss to produce good binary codes (LquantL_quant) are computed for updating of parameters of the encoder, decoder and prior (see Equation 6). (b) In the generation phase, the binary latent variables are sampled from the Quantum or classical Boltzmann Machine (BM) prior and transformed into the latent matrix through the NTN2. The output SMILES sequences are predicted by autoregressive manner. We examined the effect of the molecular generation model using the NHF. This model inputs and outputs the SMILES sequences as representations of chemical structures. We trained the model using the ChEMBL public dataset, and generated the compounds from the BM prior and the Decoder. We evaluated the validity of the generated SMILES, a fraction of SMILES strings that were successfully converted into molecular structures, as a metric for generative models. Our NHF outperformed the validity of the Gumbel-Softmax binarization (52.2% in Gumbel-Softmax and 62.0% in NHF for the fully classical computations) (Table 1). Table 1: Metrics of compounds generated by Transformer-based encoder-decoder models.00footnotetext: Validity is a fraction of the SMILES sequences that can be interpreted to the molecular graph by employing RDKit. Uniqueness is a fraction of compounds that are uniquely (not redundantly) found in the valid compounds. Sampler Binarization Validity (%) Uniqueness (%) Unique compounds / generated compounds Classical BM Gumbel-Softmax 52.20 99.94 5126 / 10000 NHF 61.95 98.09 6077 / 10000 QBM Gumbel-Softmax 71.88 95.10 6835 / 10000 NHF 96.97 51.92 5035 / 10000 2.4 Employment of the quantum annealer for sampling latent variables from the prior The major challenge in training BM is sampling from the model distribution. The Quantum Boltzmann Machine (QBM) [amin_quantum_2018] is the learning scheme of the quantum Boltzmann distribution derived from a transverse field Ising Hamiltonian using the quantum annealing processor (see Equation 1). To the effect of QBM, we trained the generative model with QBM or classical BM (see Methods). The loss of the QBM model converged within 300 epochs as the classical BM model did (Figure S3). As a consequence, the validity of compounds generated by QBM was higher than those from classical BM (97% in QBM and 73% in classical BM) (Table 1). Next, we compared the distributions of molecular properties and drug-likeness score (QED score [bickerton_quantifying_2012]) from generated compounds. Surprisingly, compounds sampled from QBM shifted their distribution of the QED score towards the higher QED side. Moreover, the proportion of drug-like compounds (QED > 0.7) was higher than those from classical BM and even the training data (Table 2). This phenomenon was specific to QBM, while samples from classical BM simply reproduced the training data (Figure 2). It should be noted herein that the afore-mentioned, molecular generation beyond the training dataset in the quality related to the drug-likeness score was achieved without any inductions of the trends, such as some restraints and stochastic conditions. Figure 2: Distributions of various molecular properties of generated compounds. Three properties are shown: drug-likeness score (QED), molecular weight (MW), lipophilicity (ALOGP), and synthetic accessibility score (SAScore). All of these properties were calculated by employing RDKit. Gray bars denote the training dataset, and orange, cyan, blue and pink lines denote the samples from the Gumbel-Softmax with classical BM, the Gumbel-Softmax with QBM, the NHF with classical BM, and the NHF with QBM, respectively. Table 2: Drug-likeness of generated compounds by Transformer-based encoder-decoder models.00footnotetext: We defined the drug-like compounds as those with the QED score of 0.7 or higher and showed the drug-like compounds as a proportion of unique compounds. Data Sampler Binarization Drug-like compounds (%) Training data 31.61 Generated samples Classical BM Gumbel-Softmax 40.81 NHF 43.15 QBM Gumbel-Softmax 52.60 NHF 66.79 Notably, under the QBM prior, the NHF-based model showed superior validity (Table 1) and drug-likeness (Figure 2, Table 2) compared to the Gumbel-Softmax-based model. Through the training steps, the coupling weights between the latent variables (spins) (see Equation 1) are optimized along with the other parameters. At the end of the training, the NHF-based model acquired denser spin-spin coupling than the Gumbel-Softmax-based model (Figure 3). Figure 3: Coupling weights of the Ising Hamiltonian trained in the generative models. (a)-(b) Heatmaps of the coupling weights between visible and visible units (Jij\J_ij\ in Equation 27; a) and between visible and hidden units (Lik\L_ik\ in Equation 27; b) of Gumbel-Softmax-based model with QPU. Note that J is the upper-triangular matrix and there are no interactions between hidden units. (c)-(d) Heatmaps of the coupling weights of NHF-based model with QPU. Subsequently, we compared relationships between transformation of chemical structures and improvement of the QED scores in the training data and samples generated by the QBM model (i.e., the latter is referred to herein as the generated data). We picked up 42 representative compounds from the training data and collected a subset of compounds in the generated data, for which substructures were matched with those of the representative compounds in the training data. This subset of the generated data is further distinguishably comprised of characteristic subsets using the following two similarity metrics: the Identically Assigned Coefficient (IAC) of matched heavy atoms and the Tanimoto Similarity (TS) of the molecular fingerprints (see Methods for details) (Figure 4). Figure 4: Relationships between structural transformation and QED improvements via QBM models. (a) The scatter plot of Identically Assigned Coefficient (IAC) and Tanimoto Similarity (TS) with respect to the subset of generated compounds, that was selected from the QBM model so as for substructures of generated compounds to be matched with those of the training data, as described in the main text (see Results and Method sections). The red and blue points denote the regions with high IAC and low TS (IAC > 0.5 and TS < 0.2), and low IAC and high TS (IAC < 0.5 and TS > 0.4), respectively (see Methods for detail). (b)-(c) Differences of the QED score between the generated and original (matched) compounds in the generated and training data, respectively. In the subset of compounds, the features are defined by “is_QED_improved”, which means that the difference of the QED score is larger than zero. The numbers in the histograms denote those of compounds that exhibit improved QED score values (red in (b) or blue in (c)) and not improved QED (gray), respectively. More specifically, in the scatter plot of the IAC and TS (Figure 4), the region with high IAC and low TS values (IAC > 0.5 and TS < 0.2) includes the following compounds: the chemical structures significantly differ from the matched, original compounds in the training data, while the molecular sizes are comparable to them. For example, a chain group of the original compounds in the training data is converted into ring structures, which may thus lead to a type of “scaffold hopping” from the original compounds. Conversely, the region with low IAC and high TS values (IAC < 0.5 and TS > 0.4) contains compounds with their scaffold conserved as that of the matched, original compounds in the training data (“scaffold preserving”), although their side chains are modified. We calculated the difference of the QED score between the generated and original (matched) compounds in the generated and training data, respectively. As a consequence, we found that the high-IAC and low-TS subset (the “scaffold hopping” subset) exhibited a higher proportion of compounds with improved QED compared to the low-IAC and high-TS subset (the “scaffold preserving” subset) (see Figure 4). The frequencies at which some substructures are found differ between generated data and training data; for example, ether substructures are found in training data at a frequency of ∼ 50%; whereas in generated data such substructures are found at a frequency of ∼ 26%. Since ether substructures exhibit some trends in terms of molecular properties (e.g., ether may improve the solubility of some compounds), we examined whether ether substructures could affect the QED-improvement distributions by comparing the QED improvement rates of compounds with ether substructures (see Figure 4) versus those compounds without ether substructures (see Figure S4). The analysis showed that compounds with ether substructures versus those compounds without ether substructures only marginally affected the QED-improvement distribution for generated data. In fact, the distribution of QED improvement was almost identical for those compounds with ether substructures as for those without, suggesting that the shift in the QED score was not caused by any substructure bias in the molecules generated by the QBM model. This feature of our quantum computational system with the QBM prior and NHF-based loss function is very useful for drug design, since scaffold hopping is still actually an intractable task. 2.5 Adaptation of the NHF to an MLP-based architecture To evaluate the scope of applicability of the NHF, we additionally performed the compound generation employing an MLP-based encoder-decoder architecture. Using a 3-layer encoder and decoder, the validity of the MLP-based model with the NHF (the NHF validity) was also superior to that of Gumbel-Softmax (the Gumbel-Softmax validity) (Table 3). Compared with the afore-mentioned, Transformer-based architecture, MLP-based architecture increased the uniqueness of generated molecules. Thus, with a modest reduction in validity, the number of unique compounds is comparable to that of the Transformer-based architecture. Note that the Gumbel-Softmax validity coupled with the QBM model was as small as 38.0%, which means that the MLP architecture does not work without employing NHF as a molecular generation model. Table 3: Metrics of compounds generated by MLP-based encoder-decoder models.00footnotetext: Definitions of validity and uniqueness are identical to those of Table 1. Sampler Binarization Validity (%) Uniqueness (%) Unique compounds / generated compounds Classical BM Gumbel-Softmax 38.50 100.00 3848 / 10000 NHF 58.30 74.90 4367 / 10000 QBM Gumbel-Softmax 39.30 100.00 3929 / 10000 NHF 59.90 86.30 5198 / 10000 3 Discussion 3.1 Role of the NHF in our objective functions (I): beyond VAE in posterior approximation We developed a novel training scheme for autoencoder-based generative models in which the latent variables are converted from continuous to discrete by the NHF for integration with quantum annealing BM, as the prior. As a consequence, the performance related to the molecular generation definitively increased with respect to the two distinct types of the autoencoder architecture (i.e., Transformer- and MLP-based autoencoders) which suggests that our NHF is a robust scheme applicable to a wide range of architectures. Compared to the loss function of DVAE (i.e., LelboL_elbo in Equation 3), the KL-divergence term is replaced with the binarization loss in the NHF (Equation 6), which thereby minimizes the gap between the binary and original (continuous) latent variable vectors. In this manner, a novel loss function is adopted in our present generative models. Although in DVAE, KL-divergence works as the regularization term, which makes the posterior distribution closer to the prior distribution (i.e., posterior approximation), balancing both reconstruction error and regularization is difficult due to repulsive effects on the entire objective. Most proposed methods for alleviating the shortcomings of the conventional VAE scheme derived from KL-divergence [higgins_beta-vae_2017, tolstikhin_wasserstein_2017, zhao_infovae_2017] require careful tuning of hyperparameters. In contrast, our presented approach is simpler: remove the KL-divergence term and establish the binarization of the latent variables, resulting in a constraint for the generative AI system (i.e., requiring the conversion of discrete and continuous latent variables). Thus, let us consider the meaning of our NHF loss function from a probabilistic perspective. We assume that the continuous vectors =fϕ()∈ℝdh=f_φ(x) ^d as outputs from the encoder fϕf_φ as well as the binary vectors ∈0,1dz∈\0,1\^d are stochastic variables. We start from a minimization of the KL-divergence between an “approximated” joint posterior qϕ(,|)q_φ(z,h|x) and a “true” posterior pθ(,|)p_θ(z,h|x) [kingma_auto-encoding_2013], of which the latter is postulated to involve the “true” distribution of data, as follows: DKL(qϕ(,|)∥pθ(,|)) D_KL (q_φ(z,h|x)\ \|\ p_θ(z,h|x) ) =qϕ(,|)[logqϕ(,|)pθ(,,)]+logpθ() =E_q_φ(z,h|x) [ q_φ(z,h|x)p_θ(x,z,h) ]+ p_θ(x) (14) logpθ()−DKL(qϕ(,|)∥pθ(,|)) p_θ(x)-D_KL (q_φ(z,h|x)\ \|\ p_θ(z,h|x) ) =qϕ(,|)[logpθ(,,)qϕ(,|)]≔ℒ =E_q_φ(z,h|x) [ p_θ(x,z,h)q_φ(z,h|x) ] (15) where ℒL is the ELBO to be maximized. Next, we postulate that the joint probability pθ(,,)p_θ(x,z,h) and the approximated posterior qϕ(,|)q_φ(z,h|x) in Equation 15 are factorized as follows: pθ(,,) p_θ(x,z,h) =pθ(|)pψ()p(|) =p_θ(x|z)p_ψ(z)p(h|z) (16) qϕ(,|) q_φ(z,h|x) =q(|)qϕ(|)=qϕ(|). =q(z|h)q_φ(h|x)=q_φ(z|x). (17) From Equation 16, the following holds: pθ(|)p_θ(x|z) is the likelihood of decoder outputs ^=softmax(gθ()) x=softmax (g_θ(z) ), which is modeled by a categorical distribution pθ(|)=Cat(|^)p_θ(x|z)=Cat(x| x). pψ()p_ψ(z) is the prior distribution of z, which is modeled by the Boltzmann Machine. p(|)p(h|z) is the prior distribution of h. From Equation 17, the approximated posterior is expressed by a single term qϕ(|)q_φ(z|x) because the encoder fϕ:↦f_φ:x is a “deterministic” function. We define p(|)p(h|z) in Equation 16 and qϕ(|)q_φ(z|x) in Equation 17 as follows: p(|) p(h|z) ≔1(2π)dexp(−λfro2‖−‖2) 1 (2π)^d (- _fro2\|h-z\|^2 ) (18) qϕ(,|) q_φ(z,h|x) ≔∏i=1dq(z(i)|h(i)),q(z(i)=1|h(i))≔1h(i)≧00h(i)<0, _i=1^dq(z^(i)|h^(i)), 14.45377ptq(z^(i)=1|h^(i)) cases1&h^(i) 0\\ 0&h^(i)<0\\ cases\ \ , (19) where z(i)∈0,1z^(i)∈\0,1\ and h(i)∈ℝh^(i) are i-th element of z and h, respectively. We maximize the expectation of the ELBO ()ℒE_D(x)L, where ()≔1N∑n=1Nδ(−n)D(x) 1N _n=1^Nδ(x-x_n) denotes the empirical distribution represented by employing the delta functions: ()ℒ _D(x)L =()qϕ(,|)[logpθ(,,)qϕ(,|)] =E_D(x)E_q_φ(z,h|x) [ p_θ(x,z,h)q_φ(z,h|x) ] =()qϕ(,|)[logpθ(|)]+qϕ(,)[logpψ()] =E_D(x)E_q_φ(z,h|x)[ p_θ(x|z)]+E_q_φ(z,h)[ p_ψ(z)] +qϕ(,)[logp(|)]−()qϕ(,|)[logqϕ(,|)]. \ +E_q_φ(z,h)[ p(h|z)]-E_D(x)E_q_φ(z,h|x)[ q_φ(z,h|x)]. (20) Notably, qϕ(,)≔1N∑n=1Nqϕ(,|n)q_φ(z,h) 1N _n=1^Nq_φ(z,h|x_n) is identical to the “aggregated posterior” proposed in previous studies [makhzani_adversarial_2015, hoffman_elbo_2016]. Although Equation 20 is the decomposition of the ELBO similar to conventional VAEs [kingma_auto-encoding_2013, rezende_stochastic_2014], we introduced the continuous value h and the binary value z as latent variables, thereby accounting for these latent variables in our formulation. The first and second terms of Equation 20 are rewritten as follows: ()qϕ(,|)[logpθ(|)] _D(x)E_q_φ(z,h|x)[ p_θ(x|z)] =1N∑n=1Nlogpθ(n|n) = 1N _n=1^N p_θ(x_n|z_n) =1N∑n=1Nlog(Cat(n|^n))=−1N∑n=1NCE(n|^n) = 1N _n=1^N (Cat(x_n| x_n) )=- 1N _n=1^NCE(x_n| x_n) (21) qϕ(,)[logpψ()] _q_φ(z,h)[ p_ψ(z)] =1N∑n=1Npψ(n). = 1N _n=1^Np_ψ(z_n). (22) The third and fourth terms of Equation 20 are rewritten by using Equations 18 and 19 as follows: qϕ(,)[logp(|)] _q_φ(z,h)[ p(h|z)] =−λfro2N1N∑n=1N‖−‖2+const. =- _fro2N 1N _n=1^N\|h-z\|^2+const. (23) ()qϕ(,|)[logqϕ(,|)] _D(x)E_q_φ(z,h|x)[ q_φ(z,h|x)] =1N∑n=1N∑∈0,1dqϕ(,|n)logqϕ(,|n) = 1N _n=1^N _z∈\0,1\^dq_φ(z,h|x_n) q_φ(z,h|x_n) =1N∑n=1N∑i=1d∑∈0,1dq(zn(i)|hn(i))logq(zn(i)|hn(i))=0. = 1N _n=1^N _i=1^d _z∈\0,1\^dq(z_n^(i)|h_n^(i)) q(z_n^(i)|h_n^(i))=0. (24) Taken together, we obtain the following equation: ()ℒ _D(x)L =−1N∑n=1NCE(n|^n)+1N∑n=1Npψ(n)−λfro2N∑n=1N‖−‖2+const. =- 1N _n=1^NCE(x_n| x_n)+ 1N _n=1^Np_ψ(z_n)- _fro2N _n=1^N\|h-z\|^2+const. =−1N∑n=1NCE(n|^n)+1N∑n=1Npψ(n)−λfro2N‖−‖F2+const. =- 1N _n=1^NCE(x_n| x_n)+ 1N _n=1^Np_ψ(z_n)- _fro2N\|Z-H\|_F^2+const. (25) As a result, the negation of Equation 25 fundamentally matches the NHF loss (Equation 6), although the orthogonalization term is absent. Still, if we consider the weight matrices of the MLP layers in the encoder =lW=\W_l\ as stochastic variables, the appropriate prior distribution term p()p(W) is multiplied with the joint prior (Equation 16), providing the orthogonalization constraint [duan_bayesian_2020]. Thus, Equation 25 exactly matches our NHF loss function (Equation 6). It should be noted that the Frobenius norm term (Equation 6) corresponds with the cross-entropy between qϕ(,)q_φ(z,h) and pθ(,):qϕ(,)[logpθ(,)]p_θ(z,h):E_q_φ(z,h)[ p_θ(z,h)] (Equation 25). In other words, the minimization of the Frobenius norm term forces the aggregated posterior to match the prior, and vice versa. Since conventional VAE models the posterior distribution using the mean-field approximation, that is, qϕ(|)=∏n=1Nqϕ(|n)q_φ(z|x)= _n=1^Nq_φ(z|x_n) in Equation 2, the objective (in particular, the KL-divergence term) independently forces each component qϕ(|n)q_φ(z|x_n) closer to the prior pψ()p_ψ(z). Such a conventional protocol often leads to low expressiveness of the latent variables (i.e., posterior collapse) [bowman_generating_2016, serban_hierarchical_2017]. Moreover, for the implementation of the aggregated posterior, previous reports introduced other architecture types, such as adversarial networks [makhzani_adversarial_2015], which increase the complexity of the system, requiring careful control of the training procedure. In contrast, our NHF scheme requires only a natural assumption on the joint prior distribution, based on the Euclidean distance between the continuous input h and the binary code z (Equation 18), which is equivalent to the similarity preserving property in deep hashing (Equation 6). Thus, the deep hashing and aggregated posterior schemes are unified in our NHF scheme. Furthermore, our NHF scheme leads to the collection of the latent variables through the Frobenius norm in the objective function (see the following section), which corresponds to a feature driven by an empirical distribution represented by delta functions in the aggregated posterior scheme (Equation 25) [makhzani_adversarial_2015, hoffman_elbo_2016]. This characteristic of our NHF scheme enables us to use a simpler and more robust technique for deriving the aggregated posterior scheme as a combination of the posterior ∑n=1Nqϕ(|n) _n=1^Nq_φ(z|x_n) and the prior. 3.2 Role of the NHF in our objective functions (I): tractable differentiability A further advantage of our NHF is found in a crucial issue for AI systems, which is related to the differentiability (i.e., backpropagation). In our NHF scheme, the binarization loss is defined by the matrix norm for all samples in the minibatch (Equation 6 and Equation 25), which is effective for approximating the gradient through the binarization function (Equations 12 and 13). This is consistent with the previous report in which the empirical loss gains smoothness with a larger sample size, when the gradients are approximated by the Straight Through Estimator (STE) [yin_understanding_2019]. In this way, handling the collection of the latent variables not only directly makes the aggregated posterior and the prior closer but also reduces the deviation of the gradient with respect to the objective function, which is expected to guide the optimization in a more appropriate direction in the state space of the system defined by the objective function. 3.3 Role of the NHF in our objective functions (I): free from reparameterization in VAE Comparing binarization schemes reveals another advantage of the NHF. In the typical VAE, the model distribution of the posterior is limited to the specific family for applying the reparameterization. To apply the reparameterization, we search for a distribution qϕ(|)q_φ(z|x) that can be decomposed by the differentiable transformation ξϕ(,) _φ(x, ρ) and auxiliary variable ∼p() ρ p( ρ). Herein, three approaches were suggested [kingma_auto-encoding_2013] for the posterior distribution: (1) tractable and differentiable inverse cumulative distribution function (CDF); (2) computation employing a “standard” distribution and its shape parameters; and (3) their composition. Gaussian distribution is the simplest example used in almost all VAE cases. In the discrete cases, posterior distribution is modeled by a factorized Bernoulli distribution and Gumbel-Softmax (Equation 5) or the other variations using the sigmoid function for smoothing [rolfe_discrete_2016]. However, such treatments can cause the vanishing gradient. We showed that NHF models have denser connections between latent variables (Figure 3), suggesting that the NHF can pass the gradient more efficiently than smoothing-based methods. In this manner, our NHF scheme exhibits a great advantage in the training of the AI system through reducing the restriction of the distribution shapes found in VAE. We trained our models using the same dataset as a published model with the conventional VAE [gircha_hybrid_2023] in this work, and revealed that the validity of the generated compounds was improved (54% and 97% in the previous and present reports, respectively). Our NHF exhibits an appropriate affinity with QBM, which can produce the distribution distinguished from that produced by the classical BM (Figure 2). In fact, we showed that QBM models increased the population of drug-like compounds, even when compared to training datasets (Table 2). It should be noted herein that we used only the tokens representing compound’s chemical structures (SMILES) for training and did not provide any information concerning molecular properties for both training and generation processes of our generative models. Nevertheless, the QBM model successfully generated more drug-like compounds, as if the generative model would have known the relationship between the structure and properties. Why was it actualized? A potential explanation is as follows: Because QBM-driven, enhanced sampling with our NHF-based, extended objective function, which involves novel posterior and prior representations, could provide more reasonable probabilistic features, the QBM model may acquire the latent features originating from a more generalized, stochastic distribution and, thus, lead to a more appropriate distribution than the empirical distribution represented by the training dataset. Consequently, our results shed light on a new generative-model approach beyond the variational inference leveraged by sampling from a quantum annealer. 3.4 Role of the Boltzmann Machine as a prior distribution Most of the VAE models adopt a fixed distribution without learnable parameters such as the standard Gaussian distribution as the prior. The objective function of VAE tends to push the encoder to zero variance (i.e., delta function), which thus leads to “posterior holes” in the latent space, thereby resulting in the low posterior probability and high prior probability [rezende_taming_2018]. In the experiment, these holes are scattered in the wide region of the latent space rather than partitioned, and so the reconstruction from the holes generates bad samples [rosca_distribution_2018]. One possible solution to avoid them is to use a more flexible and trainable prior. Actually, in some reports the mixture of Gaussian with the trainable parameters was employed with respect to each component and coefficient of the prior [tomczak_vae_2018]. However, in most of these cases, the factorial Gaussian was applied for each component, where no correlation among the latent dimensions was postulated. In this work, we adopted Quantum BM as a prior to resolve the afore-mentioned issues, which thereby extended the following two features: 1) BM naturally incorporates the correlation of the dimension via the coupler weights Jij\J_ij\, and 2) the quantum annealer accelerates the extended sampling for the training of the correlation (Figure 3). Thus, the Quantum BM has a capability to effectively capture features of the data into its parameters and generates samples that improve the performance. 3.5 Perspectives on adaptation of quantum-driven generative models to drug discovery Computer-aided molecular design is promising for efficiently finding novel candidates appropriate for pharmaceutical drugs. Our hybrid quantum-classical generative models can generate molecules that are preferable in terms of the drug-likeness (QED score), even when compared to training datasets. Moreover, our hybrid quantum-classical generative models can also provide candidate molecules with significantly drastic structural changes (e.g., from chain to ring) in generated molecules, and, therefore, such models can be utilized in searching for novel structures via scaffold hopping, as discussed earlier (Figure 4). In this work, we demonstrated that such complicated structural modifications can be realized by using sequence-based molecular representations like SMILES. Recently, many reports about exploring representations, such as graph representations [lv_meta_2024, lv_meta-molnet_2025], that are suitable for molecular generation tasks4 have been published, indicating that with further improvements in generative AI performance, other alternative representations could be effective in the future. In actual drug discovery, the scaffold is frequently fixed and functional groups simply modified because scaffold hopping is still an intractable issue, making it difficult to realize significant improvements in molecular design. Furthermore, because few clues for constructing better scaffolds have been identified, exploring the scaffold is challenging. In contrast, the quantum-driven generative models (machine/computer) shown in this work provide us with candidates that can significantly improve various properties related to drug-likeness via substantial modifications of the training dataset’s chemical structures. Then, human beings (researchers/persons) can further modify the obtained molecular structures comprising a new training dataset that can be fed back into the quantum-driven generative models. Such a cyclical molecular-design process via computers and human beings, referred to here as computer-human cyclical molecular design, is a crucial and novel drug discovery process, thus leading to a dialogic interactive cooperation of quantum-driven generative AI models and human beings in the near future. 3.6 Concluding remarks We created a novel objective (loss) function that satisfies both differentiability and binarization for our quantum generative models. This loss function is free from both the mean-field approximation and assumptions of a type of stochastic distributions (notably, both of which are imposed in VAE), thus resulting in the improvement of outcomes derived from the VAE architecture. As a result of the analysis beyond VAE, the molecular properties of the obtained data generated by our present quantum generative models outperform even those of the training dataset. Thus, hybrid quantum and classical generative modeling is a promising technique in the near future of drug discovery field. 4 Methods 4.1 Training dataset We used the subset of ChEMBL dataset provided by [gircha_hybrid_2023] (https://zenodo.org/records/7827952). The compounds are divided into training, validation, and test sets by 8:1:1 (128,800 for training, 15,360 for validation, and 15,360 for test sets). The identifiers of the split are included in Supplementary Information. The canonical SMILES strings are obtained from the molecular structures of the dataset. 4.2 Autoencoder We employed the encoder-decoder architecture of Transformer [vaswani_attention_2017] for featurizing and recomposing the SMILES strings (Figure 1). Each token from SMILES is embedded in a dmodeld_model-dimensional vector and positional encoding was added. The embedded tensors are fed into the Transformer encoder layers. The outputs of the transformer blocks are flattened by the Neural Tensor Network layer (for details, see the Neural Tensor Network part in the Methods section) to obtain the fixed-length D-dimension vectors and projected into the dvd_v-dimension by a MLP layer with batch normalization for the latent variables, where dvd_v is the number of the visible units of the Boltzmann machine (see Section 4.4 part in Methods). In the decoder block, the binary latent variables are embedded into the D-dimension continuous vector by an MLP layer with batch normalization and transformed into m-by-dmodeld_model matrices for the key and the value of the cross-attention layer. SMILES tokens, which form the embedded tensors, are passed through both the self-attention layer with the subsequent masks and the cross-attention layer as the query. Finally, SMILES sequences are reconstructed from outputs of the decoder through the softmax layer. In the generation phase, the latent variables sampled from the prior distribution are fed into the Neural Tensor Network and decoder, and the output tokens are predicted in an auto-regressive manner. 4.3 Metrics for the evaluation of generative performance We generated the tokens from trained models and transformed them into SMILES sequences. The validity was computed by dividing the number of SMILES sequences that can be interpreted to the appropriate molecular graph by the total number of SMILES generated. We employed RDKit (https://doi.org/10.5281/zenodo.591637) for constructing molecule structures from SMILES. The uniqueness was computed as a fraction of compounds that are uniquely (not redundantly) found in the valid compounds. We used RDKit to calculate the QED score, molecular weight (MW), and lipophilicity (ALOGP); we used the RDKit Contrib directory’s SA_Score descriptor (https://github.com/rdkit/rdkit/tree/master/Contrib/SA_Score) to calculate the synthetic accessibility score (SAScore) [ertl_estimation_2009]. 4.4 Boltzmann machine as a prior Whereas priors are conventionally fixed to a standard Gaussian distribution, we parameterize the prior with Boltzmann machine (BM)-based distribution powered by the quantum annealer. The BM [ackley_learning_1985] can learn a complex multi-modal probability distribution and is thus an attractive approach for integrating quantum computing into deep generative models. The probability distribution of the BM is pψ()=exp(−βEψ())ψ,ψ=∑exp(−βEψ())Eψ()=∑i<jJijzizj+∑i=1dhizi splitp_ψ(z)&= (-β E_ψ(z) )Z_ψ, 21.68121ptZ_ψ= _z (-β E_ψ(z) )\\ E_ψ(z)&= _i<jJ_ijz_iz_j+ _i=1^dh_iz_i split (26) where ∈0,1dz∈\0,1\^d is the binary latent variables, β is the inverse temperature parameter, and ψZ_ψ is the partition function. ψ=,ψ=\J,h\ are the trainable parameters: JijJ_ij denotes the interaction coefficients between unit i and j, and hih_i denotes the bias term. The Restricted Boltzmann Machine (RBM) is a subtype of BM, widely used because of its efficient training [hinton_fast_2006]. In the RBM, a latent variable is divided into the visible units =vii=1dvv=\v_i\_i=1^d_v and hidden units =uii=1duu=\u_i\_i=1^d_u, where the interactions are limited to between visible and hidden units. We considered another BM formulation wherein the visible units form interactions, similar to the Semi-restricted Boltzmann Machine [osindero_modeling_2007]. The Hamiltonian is defined as Eψ(,)=∑i<jJijvivj+∑i=1dv∑k=1duLikviuk+∑i=1dvbivi+∑k=1duckuk.E_ψ(v,u)= _i<jJ_ijv_iv_j+ _i=1^d_v _k=1^d_uL_ikv_iu_k+ _i=1^d_vb_iv_i+ _k=1^d_uc_ku_k. (27) Since the hidden units remain independent, the conditional distribution of the hidden units can be calculated independently by p(uk|)=tanh(ukλk()),λk()=ck+∑i=1dvLikvi.p(u_k|v)=tanh (u_k _k(v) ), 21.68121pt _k(v)=c_k+ _i=1^d_vL_ikv_i. (28) We applied the distribution marginalized over the hidden units pψ()=∑(−βEψ(,))ψ,p_ψ(v)= _u (-β E_ψ(v,u) )Z_ψ, (29) to the prior distribution. Therefore, the KL-divergence term of the objective function is given by DKL(qϕ(|)∥pψ())=qϕ(|)[logqϕ(|)]−qϕ(|)[logpψ()],D_KL (q_φ(z|x)\ \|\ p_ψ(z) )=E_q_φ(z|x)[ q_φ(z|x)]-E_q_φ(z|x)[ p_ψ(z)], (30) and the gradient of the second term is calculated by the difference of energy between samples from the data (positive phase) and from the prior (negative phase) as qϕ(|)[∂logpψ()] _q_φ(z|x) [∂ p_ψ(z) ] =∼qϕ(|)[∂Fψ()]−,∼pψ(,)[∂Eψ(,)] =E_z q_φ(z|x)[∂ F_ψ(z)]-E_v,u p_ψ(v,u)[∂ E_ψ(v,u)] (31) Fψ() F_ψ(z) =∑i<jJijzizj+∑i=1dvbizi−∑k=1dulog2cosh(−λk()), = _i<jJ_ijz_iz_j+ _i=1^d_vb_iz_i- _k=1^d_u 2cosh(- _k(z)), (32) where FψF_ψ is the marginalized Hamiltonian (free energy) of visible units. This equation means that the first term is calculated by the average of the latent variables from the data and that the second term is calculated using samples from the prior distribution. We used simulated annealing with Metropolis-Hastings update to approximate sampling from BM. 4.5 Quantum Boltzmann machine and its implementation We use the D-Wave quantum annealer as a source of prior samples. When the annealing time is sufficiently long with quasistatic evolution, the quantum processing unit (QPU) generates samples that approximate a classical Boltzmann distribution [amin_searching_2015]. However, the presence of a finite transverse field introduces quantum fluctuations that can lead to deviations from the ideal Boltzmann statistics. A quantum Boltzmann machine (QBM) [amin_quantum_2018], which draws samples from the Boltzmann distribution of a transverse-field Ising Hamiltonian, can be trained analogously to its classical counterpart by optimizing a variational bound on the true log-likelihood. In this study, we introduce an RBM architecture that is restricted to the Zephyr graph (i.e., an architecture native to the D-Wave QPU), thus allowing for more efficient sampling and larger dimensions than a traditional architecture, such as bipartite and clique graphs, which are not native to the Zephyr graph. In this study, the largest possible clique and bipartite graphs that could be built on a defect-free Zephyr architecture are K2(2m−1)tK_2(2m-1)t and K2(2m−1)t,2(2m−1)tK_2(2m-1)t,2(2m-1)t, respectively, where m is the grid parameter and t is the tile parameter [boothby_zephyr_2021]. The QPU used in this study has a grid parameter of 6 and a tile parameter of 4, only allowing a clique of 88 and a K88,88K_88,88 (i.e., a bipartite graph with 88 units on the visible side and 88 units on the hidden side). To utilize the entire QPU as an RBM prior and allow for a flexible number of visible and hidden units, we use chains of qubits to achieve the desired prior dimension as well as increase the connectivity between visible units. The choice of chains and hidden nodes is defined through a heuristic node-contraction method where a sparse graph is iteratively converted to a denser graph, starting with a chain length of two and increasing the length as needed. In this method, nodes are first sorted by degree. The node with the lowest degree along with its adjacent nodes, which also have the smallest degree, are chosen for contraction. When the resulting contracted graph reaches the required dimensions, the contraction stops. The node-contraction method works with any graph and especially with D-Wave’s current and future quantum annealing architectures. One advantage of this method is that it works with QPUs with imperfect graphs and allows for a granular control over the number of visible and hidden variables. In our implementation of the Boltzmann machine, the hidden units are not connected to each other, which is compatible with the Zephyr graph as there are groups of qubits that are not connected to each other. Figure S1 (see Supporting Information) illustrates the classification of qubits within the Zephyr architecture into four distinct groups, referred to as color groups. Color groups are suitable candidates for hidden units, because a color group does not have any internal connections, whereas different color groups have external connections between each other. In the Advantage2_prototype2.6 QPU, each color group contains more than 300 qubits and all hidden units are selected exclusively from a single-color group to ensure compliance with the no-connectivity constraint among hidden units; that is, selecting hidden variables from multiple color groups would violate this constraint. If the number of hidden units were to exceed 300, couplers would need to be removed to adhere to the no-connectivity condition between hidden units; however, 300 hidden units are deemed sufficient for the scope of this study. Figure S2 (see Supporting Information) shows a Zephyr graph that we used as an RBM prior, having 128 visible and 128 hidden units. Note that we utilized the entire QPU (1215 qubits) in this prior, where chains were added to increase the connectivity of visible units. When sampling from the QPU, one may need to know the effective temperature (β) at which sampling is performed. We could estimate this temperature using the maximum log-likelihood approach based on the model parameters and rate of excitations [raymond_global_2016], which are available with a D-Wave utility (https://docs.dwavequantum.com/en/latest/ocean/api_ref_system/generated/dwave.system.temperatures.maximum_pseudolikelihood_temperature.html). In this study, we scale J and b parameters of the Hamiltonian such that the QPU produces samples with an effective temperature near 1.0. 4.6 Neural Tensor Network A compound in the training data is embedded to tensor ∈ℝC×dmodelX ^C× d_model, where C is the number of token sequences and dmodeld_model is the embedding dimension. Since the sequence length varies for each compound, the tensor size is also variable. To obtain the latent vector with fixed length, it is necessary to convert the variable length tensor to the fixed length one. The simplest way is to aggregate the vectors of each length by the pooling (max, average, etc.). However, this can cause information loss concerning the length. Another method is convolution/deconvolution for converting to the fixed length vector, wherein the range of the length is fixed in advance for determining the strides. We propose a more flexible form, which is referred to here as Neural Tensor Network (NTN), using tensor products with trainable higher-order tensors. For the end of the encoder (NTN1) (Figure 1), conversion from the variable length tensor X to fixed length vector ∈ℝDh ^D is performed as follows: bα=∑i=1dmodel∑j=1dmodel∑k=1CWijα(1)XkiXkj,hα=σ(bα),b_α= _i=1^d_model _j=1^d_model _k=1^CW_ijα^(1)X_kiX_kj, 21.68121pth_α=σ(b_α), (33) where (1)=Wijα(1)∈ℝdmodel×dmodel×DW^(1)=\W_ijα^(1)\ ^d_model× d_model× D is a trainable tensor and σ(⋅)σ(·) is an element-wise activation function. For the beginning of the decoder (NTN2) (Figure 1), conversion from the latent vector ∈ℝD ζ ^D to the tensor ∈ℝm×dmodelM ^m× d_model is performed as follows: Bli=∑α=1DWαli(2)ζα,Mli=σ(Bli),B_li= _α=1^DW_α li^(2) _α, 21.68121ptM_li=σ(B_li), (34) where (2)=Wαli(2)∈ℝD×m×dmodelW^(2)=\W_α li^(2)\ ^D× m× d_model is a trainable tensor. 4.7 Derivation of the gradient of NHF with respect to the encoder parameters Using chain rules, the gradient of the loss in Neural Hash Function LnhfL_nhf (Equation 6) with respect to ϕφ is ∂Lnhf∂ϕ ∂ L_nhf∂φ =∂Lnhf∂n∂n∂n∂n∂ϕ = ∂ L_nhf _n _n _n _n∂φ =∑n=1N(∂Lrec∂n+∂Lprior∂n+∂Lquant∂n)∂n∂n∂n∂ϕ. = _n=1^N ( ∂ L_rec _n+ ∂ L_prior _n+ ∂ L_quant _n ) _n _n _n∂φ. (35) Using vector-form, LquantL_quant can be rewritten as Lquant L_quant =λfro2N‖−‖F2+λortho2‖llT−‖F2 = _fro2N\|Z-H\|_F^2+ _ortho2\|W_lW_l^T-I\|_F^2 =λfro2N∑n=1N‖n−n‖2+λortho2‖llT−‖F2 = _fro2N _n=1^N\|z_n-h_n\|^2+ _ortho2\|W_lW_l^T-I\|_F^2 (36) and the gradient of LquantL_quant with respect to nz_n are computed to ∂Lquant∂n=λfro2(n−n). ∂ L_quant _n= _fro2(z_n-h_n). (37) When the transformation from nh_n to nz_n is non-smooth, such as the Sign function or Heaviside step function, the gradient ∂n∂n _n _n cannot be computed. We approximated this derivative by the identity function. With this approximation, the backward pass is different from the forward pass, whereas the sign of negation of the backward pass is consistent with the direction that minimizes the loss [yin_understanding_2019, hinton_neural_2012]. From Equation 35 and Equation 37, we obtained ∂Lnhf∂ϕ=∑n=1N(λfro2(n−n)+∂Lrec∂n+∂Lprior∂n)∂n∂ϕ. ∂ L_nhf∂φ= _n=1^N ( _fro2(z_n-h_n)+ ∂ L_rec _n+ ∂ L_prior _n ) _n∂φ. (38) 4.8 Relationships between transformation of chemical structures and improvements of QED score We selected the representative compounds from the training data using the following criteria: 1) Molecular Weight (MW) greater than 200 and less than 450, 2) logP greater than 0 and less than 4, 3) PSA less than 150, and 4) without alert structures (computed by RDKit). From the generated compounds, we extracted those that shared the same substructures as those of each of representative compounds. As queries of substructure search, a pool of fragments is required to consist of chemically meaningful components. We created fragments from the representative compounds using RECAP [lewell_recapretrosynthetic_1998] and BRICS [degen_art_2008] fragmentation and filtered by MW (greater than 150). Herein, for a similarity measure, the Identically Assigned Coefficient (IAC) of atoms was computed as follows IAC=HACmatchedHACgenerated+HACtrain−HACmatched splitIAC= HAC_matchedHAC_generated+HAC_train-HAC_matched\\ split (39) where HACgenerated,HACtrainHAC_generated,HAC_train, and HACmatchedHAC_matched denote the Heavy Atom Count (HAC) of the generated compound, the training data, and the matched substructure, respectively. For another similarity measure, we used the Tanimoto Similarity of the molecular fingerprint. Molecular fingerprints are generated as 2048-bit vectors using Morgan algorithm [morgan_generation_1965] with radius of 2, which is roughly equivalent to ECFP4 fingerprint [rogers_extended-connectivity_2010]. 4.9 Training hyperparameters As inputs of the autoencoder, the tokens of SMILES were embedded into 160-dimension vector. For the encoder and the decoder, 5 stacks of the transformer layers with the 1024-dimension feedforward layers, the ReLU activation functions, and the Dropout with rate of 0.1 were used. The output dimension of the Neural Tensor Network 1 (NTN1) was 1024 and the output of the NTN2 was a 128-by-160 matrix. We set the number of visible and hidden units in the Boltzmann machine as 128 and 128, respectively, and trained models with the minibatch of 512 for 300 epochs, using the Adam optimizer [kingma_adam_2014] (β1 _1 = 0.9, β2 _2 = 0.999) with a learning rate of 10-4 and a weight decay of 10-3. 5 Data Availability In the current study, we used the subset of ChEMBL dataset provided by the previous report [gircha_hybrid_2023] (https://doi.org/10.5281/zenodo.7827952). The code and the generated molecules during the current study are available from the corresponding author on reasonable request. 6 References References 7 Acknowledgements This study received no funding. 8 Author contributions H.K., M.R. Y.I., V.V.C., W.K., K.C., M.A. and M.T. designed research. H.K., Y.I., Y.H., A.S. and M.T. constructed the generative models, ran experiments and analyzed the data. M.R., M.W., V.V.C., W.K., K.C. and M.A. contributed to implementation of the quantum annealing. All the authors discussed the results and contributed to writing the manuscript. 9 Competing interests The authors have no competing interests.