Hybrid Generative AI for De Novo Design of Co-Crystals with Enhanced Tabletability

Authors: nina  gubina
andrei  dmitrenko
gleb  solovev
lyubov  yamshchikova
oleg  petrov
ivan  lebedev
nikita  serov
grigorii  kirgizov
nikolay  nikitin
vladimir  vinogradov
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-14 23:09:51

Abstract: Co-crystallization is an accessible way to control physicochemical characteristics of organic crystals, which finds many biomedical applications. In this work, we present Generative Method for Co-crystal Design (GEMCODE), a novel pipeline for automated co-crystal screening based on the hybridization of deep generative models and evolutionary optimization for broader exploration of the target chemical space. GEMCODE enables fast *de novo* co-crystal design with target tabletability profiles, which is crucial for the development of pharmaceuticals. With a series of experimental studies highlighting validation and discovery cases, we show that GEMCODE is effective even under realistic computational constraints. Furthermore, we explore the potential of language models in generating co-crystals. Finally, we present numerous previously unknown co-crystals predicted by GEMCODE and discuss its potential in accelerating drug development.

Main paper: 
 introduction: 
 The use of multi-component molecular crystals, specifically co-crystals, have become increasingly popular in various industries including energy [1] , electronics [2, 3] , optoelectronics [4, 5] , food [6] , and especially in pharmaceuticals [7] [8] [9] . Pharmaceutical co-crystals are defined as solids that are crystalline singlephase materials composed of a drug molecule and an additional pharmaceutically acceptable molecule (coformer) [10] . Co-crystals have a different crystal structure from the original components, leading to unique physicochemical properties. They are appealing because the resulting solid can exhibit better physicochemical properties compared to either of the pure molecules [11] . The formation of co-crystals has been shown to enhance characteristics such as bioavailability [12, 13] , solubility [14] [15] [16] , stability [17] [18] [19] , pharmacokinetics [20, 21] , and mechanical properties [14, 22, 23] . Plasticity is a mechanical property that is particularly important for the pharmaceutical industry. It is known that highly plastic materials tend to produce stronger tablets compared to those exhibiting elastic behavior [24] . In other words, it possesses improved tabletability, defined as the capacity of a powdered material to be transformed into a tablet of specified strength under the effect of compaction pressure [25] . Therefore, it is essential to control for tabletability as it allows direct pressing with minimal addition of excipients to form a stable compact tablet. Despite all the robustness and versatility of co-crystals, determining the combination of a coformer and parent component with the desired property modification is an extremely non-trivial task, usually addressed by experimental high-throughput screening [26, 27] . Due to the large amounts of time and effort required, such studies remain targeted, focusing on rather narrow classes of candidate compounds.Artificial intelligence (AI) methods have recently found their way into the field of chemistry [28] [29] [30] [31] [32] . Since then, the accumulated experimental data has become the basis for predictive models transforming the traditional way science works. With big data and machine learning (ML), it is now possible to consider a much larger set of candidate molecules for a given problem, rather than being satisfied with a limited number of experiments. Among the pioneering works in the co-crystal domain are the studies aimed at determining the probability of co-crystallization of a particular molecular pair [33, 34] . However, the sole fact of co-crystallization with no information about the properties of the resulting co-crystals is not enough to inform decision making for a specific use case. Accordingly, another direction of research investigated co-crystal properties with AI methods [35, 36] . Still, prediction of most properties has been possible only in the case of already known co-crystallising molecular pairs. De novo design of co-crystals with predefined properties leveraging big data to cover a large chemical space remains an actual task of great application value.Therefore, here for the first time we develop a pipeline that generates coformer candidates based on the structure of a drug molecule to form a co-crystal with predefined mechanical properties. For that, we trained several state-of-the-art generative models on a dataset of 1.75M chemical structures and then fine-tuned them on the state-of-the-art dataset of coformers. We then trained a classical ML model to predict plasticity parameters of the generated coformer candidates. We further employed evolutionary optimization leveraging the trained ML models to improve the tabletability profiles of the generated coformers. Finally, we applied a pretrained graph neural network (GNN) to rank the molecular pairs according to the probability of successful co-crystal formation. We systematically evaluated and optimized the aforementioned individual components to assemble GEMCODE, a practical solution achieving state-of-the-art performance even within computational constraints. The output of GEMCODE is a set of coformers forming a co-crystal with improved tabletability properties for a selected drug compound. Thus, the pipeline can serve as a tool for selecting the best molecular combination of an active pharmaceutical agent and a coformer delivering the desired properties of the co-crystal. In essence, this work makes the following novel contributions:• We train a transformer-based conditional variational autoencoder (T-CVAE) setting the new state of the art for the coformer generation task, and hybridize it with multi-objective evolutionary algorithm to improve the desired properties of coformers. • We develop machine learning models for the prediction of mechanical properties of cocrystals for the first time in the field. • We present GEMCODE, a generative pipeline for de novo co-crystal design with target physicochemical properties contributing to drug tabletability. • In addition, we explore the capabilities of language models in the coformer generation task.• Finally, we predict a set of molecules forming novel tabletable co-crystals with known drugs. 
related work 2.1 generative ai for molecule generation: 
 Traditionally, the process of discovering new molecules or selecting chemical structures for a particular task relies on existing experimental evidence and subjective research experience, both limiting the number and variety of possible compounds to consider. Generative models allow efficient exploration of the molecular space, which has already caused a rapid growth of molecular generative design. Recurrent neural networks [37] [38] [39] [40] , variational autoencoders [41] [42] [43] [44] , generative adversarial networks [45] [46] [47] [48] , evolutionary algorithms [49] [50] [51] [52] [53] and hybrid models using reinforcement learning techniques [54] [55] [56] [57] have been successfully applied for various problems in chemistry. In this work, we trained, evaluated and compared multiple generation approaches, such as LSTM-based GAN, transformer-based VAE and conditional VAE. The latter was inspired by a study using a conditional VAE model with an attention mechanism to generate molecules [58] . However, our approach differs significantly in that we generated a condition vector based on the predictions of the pretrained gradient-boosting model. In addition, our approach includes a fine-tuning phase on a state-of-the-art dataset of coformers.DeepMind has recently presented GNoME, an AI tool for generating previously unknown inorganic crystalline materials [59] . Other similar tools exist for inorganic compounds [60, 61] . Our work also lies in the field of solid-state chemistry, but differs in the task of generating coformers, which are small organic molecules. To our knowledge, generative approaches have not yet been applied to produce coformer structures with high co-crystallization potential with drug targets. Our work effectively addresses this problem.Research in co-crystal property prediction is targeted at determining various parameters, such as the lattice energy, density, melting temperature, crystal density, enthalpy and entropy of melting, as well as ideal mole fraction solubility of co-crystals [62] [63] [64] [65] . However, a limited number of samples is typically used in the training phase. For example, Gamidi and Rasmuson trained an artificial neural network on the data of 30 co-crystal systems for 8 different drugs [35] . Such models are likely to have very limited generalization power beyond the training data. The most recent model predicting the co-crystal density [36] used a large training set of 4144 molecular pairs covering a much wider chemical space of possible co-crystals. In this work, we predict several mechanical properties of co-crystals for the first time. We use an even larger amount of data for that (6029 samples), which makes our approach more versatile and better generalizable for different pharmaceutical applications.Large language models have recently been challenged with multiple chemistry tasks, such as property prediction, yield prediction, text-based molecular design, and others [66] . The results suggest that language models are less competitive in generative tasks requiring a deeper understanding of molecular SMILES strings, but show competitive performance in classification and ranking tasks. Another study on the applicability of language models without prior specialization in the chemistry domain found that LLMs can effectively interpret chemical structures given various representations [67] . In addition, the use of language models as agents was explored in ChemCrow [68] , which makes chemistry more accessible to researchers with less domain expertise. Following up on these pioneering works, we explore the applicability of language models to the creation of coformer molecules with desired properties, which has not yet been addressed in the past. 
data collection: 
 Large dataset of molecules. In order to train a generative model capable of suggesting reasonable chemical structures, a dataset of molecules from the ChEMBL database (available with CC BY-SA 3.0 license) was collected. From the large variety of molecular structures available in the database, ∼1.75M samples were selected using criteria based on the distributions of relevant parameters in the known coformers (Appendix C.1). Using these criteria ensures that the generative models are trained on molecules capable of forming co-crystals.Dataset of coformers. Chemical structures in the ChEMBL database are still substantially different from the structures composing co-crystals. Coformers most often have more basic chemical structures and a smaller variety of functional groups. Therefore, we used an open dataset of 6819 two-component co-crystals [33] (available with MIT license), which contains 4227 unique chemical structures of the coformers, for fine-tuning.Dataset of co-crystals mechanical properties. For the mechanical properties of co-crystals, we used the Cambridge Structural Database (CSD) [69] and a recently proposed protocol for geometric analysis of co-crystalline materials available with a CSD Python API [24] . For each of the 6819 available co-crystals, we used the API to query additional experimental data from the CSD and calculate the following binary parameters of plasticity: presence of non-overlapping Miller planes (Unobstructed planes), presence of orthogonal planes (Orthogonal planes), and presence of hydrogen bonds between the planes (H-bond bridging). Since some of the co-crystals were missing in CSD, this process yielded a total of 6029 records. This data was then used for training ML models to predict each of the three plasticity parameters. We analyzed the number of samples for each plasticity parameter in the collected dataset (Appendix C.2). In the case of orthogonal planes, we observed a dramatic difference between the two groups. When training the corresponding ML model, we accounted for this disproportion by adjusting a threshold probability for predicting a positive class.Cutting-edge generative models use string [70] [71] [72] , 2D [73] [74] [75] and 3D [76] [77] [78] molecular graphs as molecular representations. The most common way is the SMILES (Simplified molecular-input line-entry system) notation, as the other approaches have not yet shaped the field to such an extent [79] . Therefore, we used the SMILES representations to describe the composition and structure of chemical molecules with short strings. Additionally, molecular fingerprints allowed us to represent molecules in a vectorized form and compare different structures by calculating a similarity measure (Appendix C.3).We used RDKit to generate 43 molecular descriptors for each coformer with its SMILES representation. Since co-crystals consist of two coformer components, each one was described by 86 numerical features in total. Before training ML models for the prediction of mechanical properties, we applied a set of preprocessing steps. We engineered new features by aggregating the molecular features of the coformers of the same co-crystal with summation and averaging. To reduce redundancy in the feature space, we investigated the feature importances using embedded methods and the degree of linear association with target variables through correlation coefficients. After feature engineering and filtering, the datasets for the prediction of non-overlapping planes, orthogonal planes, and hydrogen bonding contained 29, 24, and 30 features, respectively.Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results. 
gemcode: generative evolution-based method for co-crystal design: 
 We present GEMCODE, a novel pipeline for generative co-crystal design with improved tabletability properties. It based on the idea of hybridization of deep generative models and combinatorial optimisation. GEMCODE consists of four key components, as depicted on Figure 1 .First, a trained and fine-tuned generative model generates SMILES representations of coformer-like chemical structures. The generated molecules are then fed into the trained ML models along with the therapeutic compounds, where the mechanical properties of co-crystals are predicted. In addition, an evolutionary algorithm is used in combination with the ML models to further improve the tabletability of the generated coformers. Finally, co-crystals with the desired properties are selected for the next step, where a pretrained graph neural network scores and ranks molecular pairs of drugs and coformers according to the probability of co-crystallization. Thus, the pipeline outputs a list of potential coformers with the desired mechanical properties of the co-crystal, ranked according to the probability of successful co-crystallization. In the following sections, we describe the individual components of the pipeline in more detail.Since the number of training examples available for prediction of mechanical properties was only 6029, we resorted to the classical machine learning algorithms. We formulated a binary classification problem for each of the mechanical properties and implemented a number of ML models as a first screen, including logistic regression, k-nearest neighbors classifier, support vector machines, decision trees, multilayer perceptron, as well as ensemble models, such as random forest and gradient boosting.We then selected the best models and optimized their hyperparameters to achieve top performance. Those pretrained models were then integrated into the coformer generation and the evolutionary optimization frameworks. To validate this solution, we used an AutoML tool to design the modeling pipeline in an automated way (details are provided in Appendix G.5.The performance of a particular deep neural network is largely determined by its architecture, as well as the strategy to learn the hidden representations [80] . In order to find the most effective solution for the coformer generation task, we implemented and systematically compared three different architectures. Our evaluation included a GAN model with recurrent neural networks for both, generator and discriminator, and two transformer-based models implementing a VAE. For more information regarding the model architectures, refer to Appendix D.4 and D.5.GAN-based methods consider molecule generation a minimax game, which consists of training a discriminator to distinguish between the real data and the samples produced by a generator (Appendix D.1). In this work, we employed an open-source GAN implementation2 using LSTM to address molecule generation as a sequence-to-sequence (S2S) problem, inspired by the work of d'Autume [81] . As an alternative, we opted for a transformer architecture [82] as a basis for a VAE, since it normally outperforms recurrent neural network architectures in S2S tasks [83] .Our objective was to produce co-crystals meeting specific tabletability requirements that translate to a set of target mechanical properties. We utilized a conditional variational autoencoder (CVAE) approach [84] to achieve this. By design, CVAE makes it possible to consider physicochemical characteristics of molecules and generate co-crystals with the desired properties (Appendix D.3 offers a more detailed description of the VAE and CVAE models). We used the aforementioned mechanical properties (unobstructed planes, orthogonal planes, and H-bonds bridging) as conditions for CVAE.In the following, we refer to this model as transformer-based CVAE (T-CVAE).Finally, we included a transformer-based VAE (T-VAE) for comparison, which does not consider any specific properties of molecules, for completeness of the analysis.To increase the quality of coformer generation, we applied a graph-based evolutionary algorithm to structures produced by the generative models. The software implementation is obtained from the self-developed GOLEM library [85] . The fitness function was designed to reinforce the mechanical characteristics of molecules based on predictions of the classification models described above:f (x) = (1 -p u (x), 1 -p o (x), p h (x)) T ,where x is an evaluated molecule of coformer, p u (x) is the probability of the positive class for unobstructed planes, p o (x) is the same probability for orthogonal planes, and p h (x) -for H-bond bridging. Therefore, minimization of the fitness function f leads to generation of coformer molecules having an improved tabletability profile. Determining the possibility of co-crystallization by molecular pairing is an important step in the co-crystal design. For this reason, many works attempted to solve this problem with AI [86, 34, 87] . Most works that are closely related to our problem do not provide code to reproduce or reuse their results [88] [89] [90] [91] . To account for the probability of co-crystallization, we applied an existing GNNbased deep learning framework, called CCGNet [33] (available with MIT license). Unlike many of the previous works, CCGNet achieves state-of-the-art performance predicting co-crystal formation while being 100% open-source and easily reproducible. With an average balanced accuracy of 98.6%, CCGNet efficiently scores and ranks coformer candidates according to the probability of co-crystal formation. Since CCGNet was originally trained on the same database of coformers, we did not perform any fine-tuning and simply integrated the model from the open GitHub repository into the pipeline.5 Experimental studies 
prediction of mechanical properties of co-crystals: 
 Implementation details. The preprocessed dataset was randomly split into train and test sets in proportion 4:1. The train set was used to optimize hyperparameters of the models with a grid search using the 10-fold cross-validation (CV). The random grid size was 500 and concerned the following parameters: learning rate, number of estimators, subsample, maximum depth of the individual estimators. The test set was used only once, to evaluate and report the performance of the optimized models. We calculated accuracy and F1 score during the CV to select the best hyperparameter set. The use of the two metrics was important given the imbalanced nature of the "Orthogonal planes" and "Unobstructed planes" target variables (Appendix C.2, Figure 3c ). To account for the disproportion, we also adjusted the threshold for the probability of the positive class by calculating precision and recall metrics. Finally, we employed SHapley Additive exPlanations (SHAP) to interpret model predictions, which is based on sensitivity analysis investigating the effect of systematic changes in feature values on the model output [92] .Results. Overall, the GB model showed the best accuracy and F1 score compared to the other models across all tasks (Figure 2 ). Despite the high accuracy for the orthogonal planes parameter, we obtained a moderate F1 score suggesting that the final model is more likely to predict the absence of the orthogonal planes. This is attributed to the disproportion in the training examples discussed earlier. Although we demonstrated a significant improvement in metrics by introducing the probability threshold (Appendix G.2) evaluating the model trained on the processed data, it was not enough to entirely resolve this issue.We optimized the hyperparameters of the Gradient Boosting (GB) model, which resulted in the performance metrics outlined in Table 10 (Appendix G.4). Furthermore, we conducted a thorough With SHAP analysis (Appendix G.3), we learned that the number of atoms among the molecular pairs forming a co-crystal is a decisive factor in the prediction of non-overlapping and orthogonal planes. In both cases, the decrease in the number of atoms in the coformer molecules significantly contributed to the presence of non-overlapping and orthogonal planes. The descriptors associated with the number of hydrogen bond donors (HBD) also had a high degree of importance. As expected, an increase in the number of HBD resulted in the hydrogen bonds forming between planes of the co-crystal.Implementation details. The performance of generative models depends on hyperparameters and random restarts [93] . A grid search was implemented to select the best hyperparameters, and multiple trainings were conducted. The generative model was focused on generating coformer-like chemical structures, so it was pretrained on the ChEMBL dataset and then fine-tuned on a dataset of coformers. The importance of fine-tuning was illustrated using t-distributed stochastic neighbor embedding (t-SNE) to visualize the datasets (Appendix D.2). To evaluate the trained models, ten sets of 10,000 molecules were generated and various indicators were calculated, including validity (defined as the percentage of chemically plausible molecules to all generated), novelty (defined as the percentage of newly generated molecules that are not contained within the training set to all generated), percentage of duplicate molecules, percentage of target coformers, and diversity. More details on these indicators can be found in the Appendix F.Results. Analyzing experimental results of coformer generation, we observed that T-VAE produced the highest percent of valid and novel molecules with by far the lowest percent of duplicated structures (Table 1 ). However, among the generated coformers, only 1.68% had the target tabletability profile, as assessed by the pretrained classification models. In contrast, when generating 10,000 candidates, T-CVAE produced 5.63% of new coformers with the required mechanical properties on average. While the diversity of target coformers3 was slightly higher for GAN, it was able to produce the intermediate 2.23% of such coformers. Therefore, we conclude that T-CVAE was the most effective approach to target coformer generation. However, the transformer architecture was also the most demanding for both, the training and the generation phases (see Appendix D.6 for more details).Ultimately, we recommend to use an ensemble of generative models whenever sufficient computational resources are available. Our findings presented in Appendix D.7 suggest that the three models produce complementary results. Collectively, GAN, T-VAE and T-CVAE generate up to 2.47 times more unique target coformers than individually.Additional experiments with language models. Inspired by the most recent applications of language models in chemistry [66] [67] [68] 94] we investigated their potential in the coformer generation task. First, we employed a reduced GPT-2 model with eight heads, four attention blocks, and 14.7M parameters. Similarly to other models, GPT-2 was pre-trained on the ChEMBL dataset and then fine-tuned on the coformers dataset (see Appendix D.9 for more details). We observed that GPT-2 produced significantly lower percent of new and valid structures per 10,000 generations (Appendix D.10). Nevertheless, the model achieved 3.32% of new molecules with the desired physicochemical properties, which is comparable to GAN and T-VAE. These results prompted us to further test a more recent and capable language model. Therefore, we trained Llama-3-8B with the low-rank adoption (LoRA) algorithm using the same ChEMBL and coformer datasets. We observed major improvements in validity, novelty and the number of duplicates among the generated molecules compared to GPT-2. Notably, Llama-3-8B produced the maximum diversity percent compared to all the tested generative approaches. However, the number of molecules with target physicochemical properties dropped to only 0.34%. Analyzing these empirical results, we conclude that language models show good potential in the coformer generation task but have to be heavily optimized to achieve competitive performance with GEMCODE. We leave this endeavour for the future work.Implementation details. The multi-objective optimization algorithm used in this work considers molecules as undirected graphs and follows the generational evolutionary scheme MOEA/D [95] .First, a population of individuals is evaluated with the fitness function. Then, MOEA/D-based selection is applied to pick individuals from the population to undergo mutation. After the variation by mutation is done, the inheritance operator is used to form the new population of individuals to proceed to the next iteration (see Appendix E.1, E.2 for more details).To choose an effective evolutionary scheme for the task we compared SPEA-2 [96] and MOEA/D (see Appendix E.5). Experiments have shown MOEA/D obtaining better results in some cases.The initial population of coformer structures (obtained with the previously described generative models) were varied by the set of mutation operators, inspired by the work of Leguy [50] . The set of mutations includes simple operations (add, delete, or replace an atom, delete or replace a bond) and more complicated, multi-step actions (delete or move a functional group, insert carbon, remove an atom if it has only two neighbors). See Appendix E.3 for more details on optimization runs.Results. To evaluate results of the evolutionary search, we compared the probabilities of coformers to possess the desired mechanical properties before ("Generated") and after ("Optimized") evolutionary optimization (Table 2 ). For that, we used the pretrained ML models to retrieve the probabilities, calculated statistics and applied the non-parametric one-sided Mann-Whitney test (see Appendix E.4, F.1 for more details). In most cases, we observed a significant increase in the median probability 4 of the target class. Notably, evolutionary optimization equalized the performance of different generative models in their ability to produce coformers with the target tabletability profile. Moreover, this process consistently yielded new coformer structures, not present in the training set or in the initial population. Rivaroxaban O=C([O-])CC(=O)[O-] YORVEJ T-VAE [98] Paracetamol C1=CC=C2C=CC=CC2=C1 LUJSIT GAN / T-VAE / T-CVAE [99] C[N+](C)(C)CC(=O)[O-] CUQKAC T-CVAE [100]In order to test the effectiveness of GEMCODE, we generated coformers for the drugs with poor ability to form a tablet by powder pressing. Among the therapeutic molecules selected for the pipeline validation were Nicorandil, Rivaroxaban and Paracetamol. For each of the listed drugs, experimentally validated molecules were found among the GEMCODE-generated coformers improving tabletability of the co-crystals (Table 3 ). More details can be found in Appendix B.1.To showcase the ability of GEMCODE to predict novel coformers with target tabletability profiles, we generated coformers for one of the therapeutic molecules, i.e., Nicorandil. GEMCODE enabled discovery of 23 unique coformer with improved mechanical properties and with the presence of functional groups as in experimentally validated tabletable co-crystals (see Table 4 in the Appendix B.2). This result demonstrates the potential of GEMCODE as an indispensable tool for accelerated drug development. Broader impact is further discussed in Appendix A.Question: Does the paper provide open access to the data and code, with sufficient instructions to faithfully reproduce the main experimental results, as described in supplemental material?Answer: [Yes]Justification: All code and data is publicly available on the GitHub repository (https: //github.com/ai-chem/GEMCODE).Guidelines:• The answer NA means that paper does not include experiments requiring code. • Providing as much information as possible in supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted. 
limitations: 
 The evidence presented above looks very promising for the practical applications of our pipeline. However, a comprehensive experimental validation involving organic synthesis of coformers and co-crystal formation followed by a tablet compression experiment is required to confirm its utility. Based on our empirical results, we anticipate the following limitations of the proposed pipeline:• The coformers molecular space may be too narrow for some applications due to the small sample size of the coformer dataset. Nevertheless, if computational power is available, it is possible to use an ensemble of generative models, which partially solves the problem by increasing the number of unique molecules generated.• Currently, the GB model is biased towards predicting the absence of orthogonal planes, leading to more false negatives in the predicted coformers. We recommend exploring an alternative set of coformers based on the other two mechanical properties only.• Low-scale screening may still result in some coformers failing to form co-crystals, particularly those optimized through evolution. Screening more coformers increases the chances of finding co-crystal pairs for a specific therapeutic agent.• While polymorphism's impact on predicting co-crystal mechanical properties is not examined here, its significance is undeniable and often understated. Despite limited reported polymorphs, their potential impact on prediction model accuracy in the co-crystal field necessitates further exploration, considering the current scarcity of polymorphism data.Most limitations of the proposed pipeline can be solved with more data available for training, which remains a major challenge for successful AI applications in co-crystallization. We are working towards collecting more data and improving its quality. Also, to date GEMCODE has been adapted mainly for pharmaceutical applications. In the future, we plan to extend GEMCODE by adding more predicted physicochemical properties and other crystal forms to be able to expand beyond the pharmaceutical field.Question: Does the paper specify all the training and test details (e.g., data splits, hyperparameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?Answer: [Yes] Justification: The experimental setting is presented in sufficient detail within the paper to understand the results. More specific information is also provided with the code.Guidelines:• The answer NA means that the paper does not include experiments.• The experimental setting should be presented in the core of the paper to a level of detail that is necessary to appreciate the results and make sense of them. • The full details can be provided either with the code, in appendix, or as supplemental material. 
conclusion: 
 In this work, we presented GEMCODE, a novel generative pipeline for de novo co-crystal design. To make it as effective as possible, we implemented the hybrid generative approach combininig positive sides of both deep learning models and combinatorial optimisation. We systematically evaluated and discussed the individual components of the pipeline achieving state-of-the-art performance in the corresponding tasks. Furthermore, we performed experiments to validate the pipeline by generating coformers for three different drugs and discovering previously unknown coformers for Nicorandil.In addition, we explored the applicability of language models in the coformer generation task and identified prospective research directions. Despite limitations associated with data availability, GEMCODE enables fast generation of unique and valid chemical structures of coformers with high probabilities of co-crystallization and target tabletability profiles. This research enhances co-crystal design for pharmaceuticals and contributes to the accelerated drug development. Thanks to data and code availability, our versatile hybrid approach might find other impactful applications in chemistry.Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?Answer: [Yes] Justification: The paper provides sufficient information on the statistical significance of the performance estimates of the GEMCODE components (see Subsection 5.3, Subsection E.5 and Subsection F.1).Guidelines:• The answer NA means that the paper does not include experiments.• The authors should answer "Yes" if the results are accompanied by error bars, confidence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper. • The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions). • The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.) • The assumptions made should be given (e.g., Normally distributed errors). • It should be clear whether the error bar is the standard deviation or the standard error of the mean. Justification: The authors thoroughly reviewed the NeurIPS Code of Ethics and ensured that all aspects of their research align with its guidelines to maintain ethical standards in their work.Guidelines:• The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.• If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 10. Broader Impacts Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided. • If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. • For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. • If this information is not available online, the authors are encouraged to reach out to the asset's creators. 
new assets: 
 Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?Answer: [Yes]Justification: New assets obtained in the article are provided in the GitHub repository (https://github.com/ai-chem/GEMCODE) with proper description and instructions for use.Guidelines:• The answer NA means that the paper does not release new assets.• Researchers should communicate the details of the dataset/code/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc. • The paper should discuss whether and how consent was obtained from people whose asset is used. • At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file. 
crowdsourcing and research with human subjects: 
 Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)?Answer: [NA] Justification: The paper does not involve crowdsourcing experiments or research with human subjects.Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor research with human subjects. • Including this information in the supplemental material is fine, but if the main contribution of the paper involves human subjects, then as much detail as possible should be included in the main paper. • According to the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector. 
institutional review board (irb) approvals or equivalent for research with human subjects: 
 Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained? 
 
 Appendix: This paper presents a new method for the generative design of organic co-crystals with the goals to advance application of machine learning to pharmaceutical co-crystal design and to accelerate and reduce cost of development of solid forms of active therapeutic molecules. Extensive experimental results and multiple case studies described in the paper provide hard evidence of the effectiveness of our approach. Therefore, we are confident that this work can have a broader impact on drug discovery and development, pharmaceutical industry in general and other related domains.While we identify the aforementioned societal impacts as strongly positive, there is a risk of malicious and unintended use, as well as inaccurate predictions affecting decision-making in the drug manufacturing process. However, we deem the potential negative impacts limited due to the complexity of regulations in the corresponding fields and the laboratory experiment being the ultimate measure of success.The validation experiment involved generation of 10000 coformers using ensemble generative models, namely GAN, T-VAE, T-CVAE. The generated candidates in combination with one of the three drugs (Nicorandil, Rivaroxaban, Paracetamol) were labeled into classes of three mechanical plasticity parameters. The criterion for getting into the final dataset was satisfaction of Unobstructed planes and H-bond bridging parameters. Experimentally validated coformers among the generated molecules were then searched for using the Index Tanimoto (IT), which is a metric of molecular structure similarity. Coformers with IT = 1 were compared with molecules from literature data and added to the Table 3 .Nicorandil. Nicorandil, a medication that dilates blood vessels, is prescribed for treating angina pectoris, a condition characterized by chest pain caused by temporary reduced blood flow to the heart muscle. Unfortunately, during the manufacturing of Nicorandil tablets, the drug can degrade chemically due to the heat produced at high compressive pressure. Generation of a set of coformers for this drug resulted in a fumaric acid molecule among them. Experimental findings have demonstrated that co-crystallizing Nicorandil with fumaric acid not only led to the success of co-crystallization but also improved tabletability properties [97] .Rivaroxaban. Rivaroxaban is an anticoagulant medication that is used to prevent blood clots Rivaroxaban is often taken orally by patients, but the drug is poorly suited for direct compression tableting. The generation of coformers for rivaroxaban using GEMCODE led to the detection of malonic acid among the resulting molecules. According to Kale et al. the formation of this co-crystal leads to its excellent plastic deformation under applied compaction pressure, resulting in successful tablet formation [98] .Paracetamol. Paracetamol is an analgesic and antipyretic from the group of anilides, but forms an unstable tablet by direct pessation. Among the coformers generated by GEMCODE was found experimentally confirmed coformers of paracetamol are naphthalene and trimethylglycine (betaine) [99, 100] . They are not only able to form a co-crystal with paracetamol, but also improves its tabletability, as demonstrated in Karki et al. and Maeno et al..For the therapeutic molecule Nicorandil, the labeled generated candidates were selected by meeting the conditions of being recognized as safe for human use and presence of carboxyl functional group, which means finding molecules forming the same synthon as experimentally confirmed coformers.The discovered coformers are presented in the for the probability of co-crystallization with Nicorandil for the candidate molecules. The ranking was performed within the set of molecules generated by each model separately.Coformer analysis. Besides carboxyl, the demonstrated chemical compounds contain functional groups such as hydroxyl and amide groups, which are characteristic of the confirmed coformers of Nicorandil from the work of Maeno et al. The generated coformers contain various structural modifications, such as changes in the length of the carbon skeleton, addition and partial substitution of functional groups, the appearance of multiple bonds and benzene rings. It is important to understand that GEMCODE is focused on the search for the best candidates from the point of view of crystallography and does not address the deep issues of interaction of the created structures with the human body. Therefore, with the help of expert evaluation, we selected molecules that should normally be safe for use in pharmaceuticals.Meanwhile, many of the discovered coformer structures not only already meet all three mechanical parameters consistent with improved tabletability properties, but also have a high probability of successful co-crystalization with Nicorandil. In this way, the process of discovering new coformers using GEMCODE can be described as a smart approach to selecting candidate molecules for propertycontrolled co-crystallization.The ChEMBL database contains information on more than 2.4 million drug-like chemical compounds.For training generative models, we needed a large number of molecular structures that would have similar properties to the known coformers. Therefore, 1.75 million samples were selected from the molecular structures of the ChEMBL database according to the following criteria:• Structural type: molecule.• Class: small molecules.• Molecular weight of each component <600 Da.• Number of hydrogen bond donors (HBD) less than 3 and hydrogen bond acceptors (HBA) less than 8.• Number of rotatable bonds up to 9.• Polar surface area up to 138 nm.• Number of heavy atoms in molecular structure up to 39.Mechanical properties of the co-crystals determine their viscoelastic nature. The presence of unobstructed planes and additional slip planes orthogonal to the stacked layers lead to the improved plasticity [101] . Also, there exists evidence that the lack of hydrogen bonding between the layers has a positive effect on the plasticity of a crystal [102, 103] . Therefore, an "ideal" co-crystal in terms of plasticity should have non-overlapping slip planes, additional orthogonal planes and no hydrogen bonding between the planes (Figure 3a ).The compaction properties of many pharmaceutical powders depend on their viscoelastic nature. The closer the material to a perfectly plastic body, the larger the bonding area after compaction of the powder and the denser (less porous) the compressed tablet is (Figure 3b ). Therefore, accurate prediction of the plasticity parameters is essential for data-driven co-crystal design.Traditionally, molecules are represented as structural diagrams with bonds and atoms, but such representations are not well suited for efficient computation. Alternatively, molecules can be represented with SMILES and molecular fingerprints, which have been extensively used for various applications, including the generative models [104] . SMILES notation is often used to describe the composition and structure of a chemical molecule by means of short strings (Figure 4a ). Whereas molecular fingerprints is a way of representing molecules in the vectorized form (Figure 4b ). Therefore, molecular fingerprints enable comparing different structures by calculating similarity measures. GANs typically consist of two neural networks, a generator and a discriminator, playing an adversarial game against each other while learning the data distribution p * (x). The generator network receives a random input signal and generates data distribution p θ (x), while the discriminator network D ϕ (x) evaluates the generated data and tries to distinguish it from the real training examples [105] . In the original formulation, both networks are improved by competing with each other following a min-max optimisation procedure: [106] :min θ max ϕ E p * (x) [logD ϕ (x)] + E p θ (x) [log(1 -D ϕ (x))].E p θ (x) [-log(D ϕ (x))].Since 2014, GANs have been successfully used for numerous applications, including modeling of astronomical phenomena [107] , experiments in particle and high-energy physics [108] , medical imaging [109] , and molecule generation [45, 46] . The GAN takes SMILES representations of molecular structures as input. In the training process, the generator network creates molecular representations from the Gaussian noise and the discriminator network tries to differentiate those from the tokenized SMILES of the real chemical compounds. As a result, the generator learns to output new molecular structures similar to those in the training set.The GAN trained on the ChEMBL dataset with batch size of 512 and learning rate of 0.001 consistently produced molecules with validity >75 % after 25,000 training steps. After 30,000 steps, this model was fine-tuned on the coformer dataset with a smaller batch size of 256 for additional 1,000 steps (Figure 5a ). The t-SNE analysis reveals that the molecular space of coformers is considerably more constrained compared to that of ChEMBL (Figure 5b ). Therefore, fine-tuning was critical to shift chemical compound generation towards the molecular space of coformers. The final model was able to produce >95% of valid and >86% of unique chemical structures molecules in the test generation of 1000 molecules at 5 times repetition. Variational Autoencoders (VAEs) consist of two deep neural networks, namely, an encoder and a decoder. The encoder network takes an input feature vector and converts it into a fixed-dimensional vector, while the decoder network converts fixed-dimensional vector back to the original input feature vector. The primary objective of autoencoder is to learn an identity function, and the fixeddimensional vector is referred to as the latent vector z. This latent vector z serves as an information bottleneck, meaning that it is designed to capture only the most statistically salient information in the data. In VAEs, the latent vectors z are sampled from a normal distribution N (0, I), where I is the identity matrix. To train VAEs, the loss function, which needs to be optimized for input data X and latent vector z, can be formulated as follows:E[logP θ (X|z)] -D KL [Q θ ′ (z|X)||N (z)],where D KL is the Kullback-Leibler divergence, which measures the difference between two probability distributions, Q and N ; E is the mathematical expectation; P and Q are probability distributions.The probability distributions P θ (X|z, c) and Q θ ′ (z|X) are learned by deep neural networks called the decoder and encoder, respectively. These networks have learnable parameters θ and θ ′ . The first term of the loss function is the reconstruction error for the input data X. In contrast, the second term measures the similarity between the probability distribution of the latent space and the target probability distribution, N (z), which is N (0, I).When using VAEs, it is difficult to control the specific properties of the generated data. Also, since the latent vector is sampled from a unimodal Gaussian distribution, the generated objects tend to be very similar to the training examples. This is not an efficient way to generate new molecules.Conditional variational autoencoders (CVAEs) were developed to address these challenges. CVAEs can learn multimodal probability distributions by adding a condition vector as an additional input during the generation process. The objective function of a CVAE with condition vector c (passed as an input to the encoder and the decoder) is given by:E[logP θ (X|z, c)] -D KL [Q θ ′ (z|X, c)||N (z)].It is common knowledge that generating long sequences can be a challenging task for recurrent neural networks. Therefore, we considered transformers, as a more modern and effective architecture for the task. In our approach, we apply the Pre-Layer Normalization Transformer [110] , a modification of the original Post-Layer Normalisation Transformer. Similarly to VAE, it is composed of two neural networks, an encoder and a decoder, but with the attention mechanism. Such architectures are known to suffer from a posterior collapse [111] . To overcome this, we used Kullback-Leibler divergence annealing (KLA) [112] . Ultimately, the loss function of the T-CVAE is given by:E[logG θ ′ (X t |z, X dec , c)] -k w D KL (Q θ (z|X enc , c)||p(z|c)),where D KL is the KL divergence; E is the mathematical expectation; Q θ is a parameterized encoder function; Q θ ′ is a parameterized decoder function (generator); p(z|c) is a conditional Gaussian prior. Here, θ , θ ′ , X enc , X dec , z , X t , c , k w are the parameter set of the encoder, the parameter set of the decoder, the input of the encoder, the input of the decoder, the latent variables, the reconstruction target, the conditions, and the weight for KLA, respectively. This objective function was inspired by the work of Kim [58] .Our proposed architectures of T-VAE and T-CVAE are shown in Figure 6 . Since T-CVAE is an upgrade of T-VAE, the schematics are virtually identical. They differ only by the presence of a block responsible for concatenating an additional condition vector with physicochemical properties with the latent space vector in yellow and the molecule token vector at the input to the encoder. Thus, the architecture represents a language model of a transformer whose encoder encodes information about molecules in a 128-dimensional conditional Gaussian latent space. In case of T-CVAE, a vector of conditions consisting of physicochemical properties predicted by gradient boosting model is attached to the latent space. Based on this vector, the transformer decoder learns to generate coformer candidates. The molecule tokens are embedded into 512 dimensions, the same as the model dimension. The encoder and decoder in the proposed transformer architecture consist of 6 layers of 8 heads each. To put this into perspective of our study, T-CVAE required about 45 minutes on an NVIDIA RTX A6000 graphics card to generate 10,000 molecules. GAN managed to do the same in 1.88 seconds. A similar evaluation on a more common NVIDIA GeForce RTX 2070 resulted in 3.73 hours and 3.37 seconds, respectively. Arguably, this makes T-CVAE practically infeasible for many users. Taking this into account, we recommend a pragmatic choice to keep GAN as the default generative model in GEMCODE. However, a combination of generative models is required to achieve broader exploration of the target chemical space.Testing whether the generative models produce the same molecules was a fascinating experiment. In our study, we performed ten generations of 10,000 molecules each, 100k molecules for each model in total. Then, we filtered molecules by discarding duplicates, chemically invalid molecules, not new molecules, and molecules that do not satisfy the required physicochemical properties and SA. After filtering, we got the following result: out of 100k molecules, GAN generated 1639, T-CVAE -2452, and T-VAE -1407 molecules. Among all of those, only 202 molecules were found to be common. A more detailed evaluation of the intersection of the generated molecules can be seen in the Figure 7 . Thus, in this experiment, each model was able to generate unique molecules that the other models did not produce. Specifically, GAN generated 1051, T-CVAE -1733, and T-VAE -698, making the total of 3482 new unique coformer molecules with the properties of interest. Therefore, taken together, the models can generate 1.42 times more molecules than T-CVAE, 2.47 times more than T-VAE, and 2.21 times more than GAN. According to these empirical results, if necessary and under certain research conditions, it may be relevant to use an ensemble of generative models to increase the total number of generations of diverse molecules. Of course, this approach is inherently resource-intensive.Table 1 shows the values of target molecules in percent. These percentages can be interpreted as the probability with which a generative model can create a new molecule with the properties of interest in a single generation iteration.In order to avoid misunderstandings, we also present below a numerical comparison with the target molecules in the training dataset.To interpret the efficiency of generative models, we present a quantitative comparison of target molecules in the training dataset. A dataset with examples of existing pairs of co-crystals, consisting of 4200 molecules, was used for additional training for the task of generating co-crystals. This dataset contains 355 molecules, which correspond to our selection conditions for the formation of a co-crystal with Theophylline. As can be seen from the description above and Figure 7 , even the least efficient model was able to synthesize 1639 new potential molecules during the experiment. Thus, it can be seen that the use of the generative models developed by us to search for a coformer is very promising. Also, it is worth noting that if we repeat the experiment and generate another 100,000 molecules, we will be able to obtain a significant number of more potential molecules. Since we have not conducted this additional study, and we determine exactly what the limit of model in generating unique molecules is, we cannot calculate any relative to compare model performance with respect to the training data. For this reason, we propose values.order to further analyze the novelty of the generated molecules, we performed additional experito calculate the Tanimoto similarities for the generated molecules. For this purpose, we plotted histograms illustrating the distribution of maximum IT between the generated coformers and the coformers the training dataset (Figure 8a ). Notably, the distribution is mostly centered on IT values between 0.5 and 0.61 for all generative models. This observation strongly supports the claim that the generated molecules exhibit substantial novelty. The Tanimoto similarity distributions between all molecules generated by GAN, VAE, and CVAE were also analyzed (Figure 8b ). For each model, the average Tanimoto similarity ranges from 0.70 to 0.75. On the one hand, this indicates a sufficient diversity of molecules. On the other hand, the relatively high average similarity was expected, since all generated coformers refer to the formation of a co-crystal with the same drug. This fact agrees well with the observation that the CVAE distribution is skewed towards 1 due to the "condition" architecture block that enhances the drug-specific targeting properties of the coformers.GPT-2 was pretrained on the ChEMBL dataset with a batch size of 128 for five epochs and then fine-tuned on the coformers dataset with a batch size of 256 for 15 epochs. We used the top-p (nucleus) sampling method with p = 0.95, as recommended in the study by Holtzman [113] for balanced quality, diversity, and coherence in the model outputs. We evaluated generations of GPT-2 and Llama-3-8B models and compared them to the best scores produced by the generative models in GEMCODE (Table 6 ). While GPT-2 produced a higher percent of target coformers compared to T-VAE (2.14% versus 1.68%), the ensemble of generative models in GEMCODE resulted in at least two times more molecules with the desired physicochemical properties (5.63% with T-CVAE only). Llama-3-8B was able to generate only 0.34%, which is not attractive given vast amount of resources required to train and use the model for generation.Meanwhile, Llama-3-8B clearly outperformed GPT-2 in terms of validity, novelty and the number of duplicates in molecule generations. It also showed almost on-par performance with GEMCODE, e.g., the percent of valid molecules (98.3%) was comparable to T-VAE (99.7%) and T-CVAE (98.4%).One important advantage of Llama-3-8B over GEMCODE is its ability to generate more diversity in the target molecular space. However, this advantage is diminished by the overall low percent of target coformers generated.To find molecules with a higher probability of exhibiting desired characteristics, we used an evolutionary algorithm based on a self-developed GOLEM framework [85] 5 . The evolutionary algorithm operates on graph representation of molecules. The problem being solved is the minimization of the multi-objective F in the discrete space of structural graphs M ð under a set of constraints C. The task is to find an optimal structure M * g = ⟨V, E⟩.M * g = argmin M ð F, where M ð = {M g | C(M g )}. F = (1 -p u (x), 1 -p o (x), p h (x)) T ,where x is an evaluated molecule of coformer, p u (x) is the probability of the positive class for unobstructed planes, p o (x) is the same probability for orthogonal planes, and p h (x) -for H-bond bridging. Therefore, minimization of the fitness function F leads to generation of coformer molecules having an improved tabletability profile.The general scheme of the evolutionary algorithm is presented in the Figure 9 . At first, individuals from an initial population are selected for mutation. At the mutation stage, the individuals are modified using a set of mutations described earlier. To control the process, Mutation operator refers to Change Advisor that determines possible actions. Since the algorithm implements a generational evolution scheme, Inheritance produces a new population using all the individuals obtained through mutation. Ellitism operator replaces the four worst individuals in the new population by the best ones found so far. The cycle is repeated until any stopping criteria are satisfied (time limit or maximal number of iterations). For each generative model, 10 independent runs of optimization were performed. For each run, we used a random sample from all unique coformers generated by the model as the initial population. The sample size was equal to the mean number of target coformers found in the 10,000 generated molecules. Analyzing the dataset of already known coformers [33] , we estimated the maximal number of heavy atoms to be 50 and the available elements to be C, N, O, F, P, S, Cl, Br, and I. Those estimates were used to configure the evolution process. Population sizes were set to 200, number of iterations to 200 and timeout to 60 minutes.Evolutionary algorithms tend to produce redundantly complicated structures due to overfitting [115] . To avoid unrealistic molecules, synthetic accessibility score (SA) [116] was calculated for all molecules obtained with evolutionary optimization. Only coformers with SA ≤ 3 were selected for further consideration.The results of the evolutionary optimization application were the most prominent for H-bond bridging and are presented in the Figure 10 . The genetic algorithms are considered to be a strong basis for drug design tasks [53, 117] , so we compare the GEMCODE with genetic baselines. Known baselines from Guacamol [118] can be used to optimise any molecule in SMILES notation with a given goal. However, unlike the Guacamol tasks, the co-crystal design task is multi-objective, whereas the algorithms from Guacamol_baselines6 (e.g. the well-known GraphGA) are focused on single-objective tasks.We have developed the multi-objective modification of GraphGA with Pareto dominance-based fitness, which can be used for the co-crystal design tasks. It started from a random subset of co-crystals (with the same population size and number of iterations as used in GEMOL). GraphGA showed inferior results according to the mechanical property values obtained from GEMOL (see Figire 12) . We can see that for unobstructed planes the best probability average over all runs is 0.95 for GEMCODE against 0.938 for GraphGA, for orthogonal planes -0.72 against 0.63 and for h-bonds bridging -0.18 against 0.12. Also, the convergence is quite unstable (we think it is caused by insufficiently successful selection procedure), so the hybrid approach implemented in GEMCODE is better in all cases.Novelty is a ratio that shows how many valid molecules without duplicates among the generated ones are novel (i.e., these molecules were not contained within the training dataset and were produced by the generative model). Therefore, Novelty is defined as follows:Novelty = N V • 100 [%].After we obtain a set of valid molecules without duplicates that are not contained within the training dataset, we predict the physicochemical properties using the pretrained gradient boosting model. We then select molecules by the required mechanical properties: Unobstructed planes = 1, Orthogonal planes = 1, H-bonds bridging = 0. For those, we also evaluate the synthetic accessibility (SA) score using rdkit.Contrib.SA_score. As mentioned before, we define the target tabletability profile in the generation process by the required mechanical properties and SA ≤ 3. Thus, the percentage of target molecules is calculated as follows:Target coformers = S G • 100 [%].F.6 Synthetic Accessibility Score (SA)SA is calculated using rdkit.Contrib.SA_Score.We used Diversity to assess the diversity of the generated molecules. In order to estimate this indicator, we need to calculate Tanimoto-similarity (T s ) and Tanimoto-distance (T d ). Thus, to estimate T s , consider two molecules, a and b, with Morgan fingerprints m a and m b , respectively. The number of common fingerprints between the two molecules is represented by |m a ∩ m b |, and the total number of fingerprints is represented by |m a ∪ m b |. Then, Tanimoto-similarity is defined by:T s = |m a ∩ m b | |m a ∪ m b | .Then, Tanimoto-distance and Diversity are related as follows:Diversity = T d (a, b) = 1 -T s .When it comes to dataset splitting, the choice of technique has a significant impact on model performance. For the field of co-crystals, there is still no established splitting strategy. The potential for more stratified splitting approaches in drug design, such as molecular scaffolds, exists, but adapting this method to co-crystal dataset presents a number of challenges. Each sample in the co-crystal dataset consists of two coformer molecules with different scaffolds and a large structural diversity, unlike drug design applications that deal with specific classes of compounds and typically allow the identification of fewer scaffolds. For out-of-distribution generalisation analysis, along with random splitting of the dataset, we separated the training and test samples based on Tanimoto similarity, maximising dissimilarity between molecules of different subsets. The results of the experiment are summarised in Table 9 . Since random splitting has already been used in the co-crystal domain within other works, it is practical and showed the best performance the random splitting was preferred. Due to the significant imbalance of data related to the orthogonal planes, we decided to change the threshold of probability of assigning a sample to a class. For this purpose we investigated the change of precision and recall depending on the value of the introduced threshold (Figure 13 ).Figure 13 : Change of metrics depending on the set threshhold.The threshold was set under the condition of equality (intersection of lines) of precision and recall metrics, as it represents the optimal point for balancing the number of false positives and false negatives.In order to increase the transparency and reliability of the ML model's decisions and results, we used the SHAP method, which allows us to interpret the output of the predictive model. But also, it enables domain-specific hypothesis generation while contributing to the explainability of the predictive model, which is a huge benefit for potential applications. The Figure 14 shows the features that were used in the training process, ranked in order of importance for the final prediction. In this case, the SHAP values to the left of the center vertical line are negative-class (0) and to the right are positive-class (1). Also, red dots indicate a higher feature value and blue dots indicate a lower feature value.One of the interesting applications of machine learning models in the field of co-crystal design is the prediction of different physicochemical properties. Table 10 presents a comparison of known models that predict properties such as crystal density, entropy and enthalpy of melting, melting temperature, ideal solubility, and lattice energy. It is important to note that the mechanistic properties of co-crystals have not been predicted before, so the metrics obtained in this paper can be considered state-of-the-art. Justification: The limitations of the proposed approach are detailed in the main text in the Section 6. Guidelines:• The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper. • The authors are encouraged to create a separate "Limitations" section in their paper.• The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be. • The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated. • The authors should reflect on the factors that influence the performance of the approach.For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon. • The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size. • If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness. • While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.https://github.com/ai-chem/GEMCODE 38th Conference on Neural Information Processing Systems (NeurIPS2024).https://github.com/urchade/molgenA numerical comparison of the target molecules generated by the models and those available in the training dataset is given in the Appendix D.7Median probability score can be best described as the median probability of assigning coformers to a positive class for each of the mechanical properties. In other words, median probability score gives an idea about the central tendency of the model's confidence in predicting a particular mechanical property.https://github.com/aimclub/GOLEMhttps://github.com/BenevolentAI/guacamol_baselineshttps://github.com/autogluon/autogluon, version 1.1.0 


Decision: This paper proposes GEMCODE: a pipeline to design novel coformers that could form tabletable co-crystals with a given organic molecule. GEMCODE is comprised of a sequence of approaches, including generative and predictive models, as well as evolutionary algorithms. The authors conduct an empirical coformer design study using known drugs, showing that GEMCODE can both rediscover known coformers and also propose promising new ones.

The paper sparked significant debate among the reviewers, both during the author rebuttal period, as well as later during the inter-reviewer discussion. Notably, one reviewer consistently championed the work (rating of “Strong Accept”), and one leaned slightly towards rejection; the other two reviewers were undecided (shifting their ratings up during rebuttal but then down during reviewer discussion) but ended up leaning slightly towards acceptance.

Main positive points raised by the reviewers were as follows:

-	This paper is one of the first efforts in its domain (subdomain of the broader “AI for Science”), and likely would be the first one accepted at a conference such as NeurIPS. Its acceptance could generate more interest in applying ML to an impactful scientific problem.
-	The paper is well-written, thorough, and makes extensive use of existing approaches/datasets where appropriate (it does not strive for novelty for novelty’s sake).
-	Experimental results, albeit purely in silico, appear promising.

Some of the issues highlighted by reviewers were addressed by authors during the rebuttal, and hence are not further discussed here. Main complaints that were not resolved are as follows (each accompanied by a counterpoint raised in discussion):

-	To confirm GEMCODE is fully practical, it would be great for the work to (a) predict/consume crystal structure instead of SMILES alone to account for polymorphism, and (b) include wetlab experiments to verify the results. While everyone generally agreed this would be ideal, the positive reviewer argued the work is already useful as a starting point for future research. However, reviewers noted that, in the presence of these limitations, the claims made in this work may appear overconfident.
-	ML approaches used in this paper are generally not novel (mostly simple adaptations of existing ideas), and the work appears more valuable from the domain side, raising questions as to whether it’s appropriate for NeurIPS. While I agree with the lack of novelty, my own comment on this would be to keep in mind that “AI for Science” is a strong theme at NeurIPS-like conferences, and while most accepted papers in this domain do have significantly more novelty than GEMCODE, they unfortunately are often of limited practical use, precisely because NeurIPS sometimes tends to select for novelty over practicality. If anything, accepting a more practical and less methodologically novel work could correct for this bias. Future works on co-crystal design would have something to compare to, and likely also bring more ML novelty.

Summing up, this paper generated a good deal of healthy disagreement between the reviewers. I think the arguments from both sides are perfectly reasonable and could see this paper falling on either side of the acceptance bar. **As the paper pioneers a new subdomain of “AI for Science” that could be of interest to the NeurIPS community, and it does so in a practical way, I am willing to lean into its positive aspects, and recommend the acceptance of the paper.**
That being said, the complaints listed above remain, and I would **urge the authors to make a serious attempt at integrating all the review feedback** into their work. In particular, please make the limitations of GEMCODE (e.g. not dealing with polymorphism, lack of experimental validation, etc) clearer, and discuss the extra references and points of view provided by the reviewers.

Reviews: review: paper_summary: paper_summary: The authors investigate an interesting chemical problem of generating coformers given an organic molecule such that they would form co-crystals with desirable chemical properties. The authors use GAN/VAE-based methods to generate SMILES of potential coformers, which are then improved by evolutionary optimization, before finally predicting the probability of co-crystallization via a GNN. Traditional statistical models appear to be able to predict certain mechanical properties of the co-crystals given the pairs of SMILES. The GAN/VAE models are able to generate valid molecules, where ~40-80% are predicted to crystallize.
strength_weakness: strength_weakness: 1. As with any ML for science work, the devil is in the details of the application. There are several big problems I see:
- Given that co-crystals are known to also have a variety of polymorphs, predicting mechanical properties based on SMILES without predicting/knowing the exact crystal structure inherently is not a sound approach (which polymorph's property is predicted?). 
- While I agree the overall crystal contacts/packing dictate mechanical properties (ignoring defects), I do not see sufficient justification as to why predicting unobstructed planes/orthogonal planes/H-bonds bridging is a strong proxy for plasticity of co-crystals, let alone tabletability (there are just so many different types of possible interactions and/or indicators). The citations there are at best weak and do not support the claims. 
- The definition of validity/novelty is very relaxed. Validity that is 'this molecule has the right valence' does not rule out thermodynamically infeasible molecules (at least, QED/SAScore should both be reported; otherwise, by definition, you can use SELFIES to get 100% valid molecules with every method). Novelty should not be 'how many molecules are not exactly the same as the training set' but rather a histogram of the closest Tanimoto similarity of the generated molecules to the training set (because a molecule can differ by 1 atom and this metric would still count it as novel). Duplicates are better represented by diversity (histogram of Tanimoto similarity between the generated molecules). Sec 5.4 and Table 4 show very similar/simple aliphatic carboxylic acids. This worries me that the molecules produced lack diversity and are not novel. 
- It is very hard to validate the results without wet lab results. The predicted probability of co-crystallization and predicted mechanical properties are unfortunately not verifiable outside of models trained here unless structures are known. While I fully agree that experimental validation should not appear in NeurIPS, I fear that 'generating co-crystals with good tableability' cannot be claimed unless there are wet lab results. 
2. The ML approaches taken here are not novel, and to an extent, questionable. If the authors developed a crystallization/property predictor (which produces a combination of scores), I feel the typical optimization approaches (e.g. reinvent, or at least the suite of software in Guacamol) should be at least used as baselines, and otherwise it is hard to justify the usefulness of VAE/GAN; the comparisons against GPT-2 models are much less relevant.

Note: citations on the chemistry side can be significantly improved (e.g. all of the citations 1-6 are hardly relevant, e.g. there are plenty of impactful reviews for charge transfer co-crystals and their applications). I also disagree with the current screening for co-crystals 'focus on rather narrow classes of candidate compounds' (especially when the demonstrated results in the paper are all aliphatic carboxylic acids).
questions: questions: 1. I note that in Sec 5.1, it says the data were randomly split. As we know, co-crystal molecules typically have limited diversity, can the authors elaborate how the predictions generalize out of distribution (e.g. by scaffold splitting)? 
2. Can the author produce novelty/diversity of generated samples as distribution of Tanimoto similarities? 
3. Could you provide CCDC access code (as is tradition) for Table 3?
limitations: limitations: yes
paper_summary: paper_summary: This paper presents GEMCODE, a novel pipeline for generating co-crystal designs with enhanced tabletability properties for pharmaceutical applications. The authors combine deep generative models, evolutionary optimization, and machine learning to create and evaluate potential co-former molecules. They train models to predict mechanical properties of co-crystals and generate co-former candidates using various approaches, including GAN, transformer-based VAE and CVAE architectures. The generated candidates are then optimized using evolutionary algorithms and ranked based on co-crystallization probability. The authors report that their T-CVAE model produced the highest percentage of co-formers with desired tabletability profiles. They validate their approach by generating experimentally confirmed co-formers for drugs like Nicorandil, Rivaroxaban, and Paracetamol. While the results appear promising, the authors acknowledge limitations such as potential bias in property predictions and the need for more comprehensive experimental validation. They also explore the use of language models for co-former generation, finding potential but noting the need for further optimization to achieve competitive performance.
strength_weakness: strength_weakness: I recommend adding these references in the introduction and related works, as they are relevant to the manuscript [1, 2, 3]. In table 1 and 2, arrows indicating the direction of optimization for the metrics would be helpful. 

While the authors provide some validation using known co-formers, there is a lack of comprehensive experimental testing of the novel co-formers generated by GEMCODE. Synthesis and physical testing of predicted co-crystals would significantly strengthen the claims about the pipeline's effectiveness; or at least some theoretical ab initio results on the validation  More details on the validation experiments are needed, in particular how the therapeutic molecules were selected. 

[1] Dollar, Orion, et al. "Attention-based generative models for de novo molecular design." Chemical Science 12.24 (2021): 8362-8372.
[2] Jensen, Jan H. "A graph-based genetic algorithm and generative model/Monte Carlo tree search for the exploration of chemical space." Chemical science 10.12 (2019): 3567-3572
[3] Tripp, Austin, and José Miguel Hernández-Lobato. "Genetic algorithms are strong baselines for molecule generation." arXiv preprint arXiv:2310.09267 (2023).
questions: questions: Can the authors confirm if the validation experimentation generated co-formers are not present in the training set (the pre-training CHEMBL set or the fine-tuning co-former set)? Are the therapeutics chosen in the validation experiments also not present in any of the datasets (the co-former set in particular)? 

The validation experiment needs to be expanded, as the results here would determine the utility of the entire GEMCODE pipeline. While the individual parts of the pipieline work, the validation experiments show that the co-crystal designs are real. Are there any fitness functions or ab initio results that can demonstrate the effectiveness of GEMCODE in generating co-crystals with novel therapeutics?

The state-of-the-art molecular generation algorithms [1], particularly for drug design tasks, are evolutionary. Rather than using language models, it seems more likely that a genetic algorithm would generate better structures than the GAN/VAEs/LLMs used here. 

[1] Tripp, Austin, and José Miguel Hernández-Lobato. "Genetic algorithms are strong baselines for molecule generation." arXiv preprint arXiv:2310.09267 (2023).
limitations: limitations: The authors acknowledge a bias in predicting the absence of orthogonal planes, which could limit the pipeline's effectiveness. More work is needed to address this imbalance in the training data or model architecture. Furthermore, while tabletability is important, the narrow focus on this property overlooks other crucial aspects of pharmaceutical co-crystals such as solubility and bioavailability. This limits the utility of the pipeline for drug development. The authors acknowledge these issues in the limitations and conclusions.
paper_summary: paper_summary: This works presents a generative framework for co-crystal design which uses deep learning and evolutionary algorithms for optimization. The GEMCODE pipeline can be used to select optimal molecular pair combinations: an active pharmaceutical, and a coformer to control for the desired co-crystal properties. The work is an application of generative AI (specifically, GANs and VAEs) to an interesting and important problem in the pharmaceutical industry. The paper has great figures and is well-written, and the code appears well-documented. Nevertheless, the paper could improve from some slight re-structuring and/or re-writing to improve the clarity of the paper, as it was tough to follow in places (there are many molecular representations and model architectures used, many target properties, these could be better organized so that it is easier to follow). But overall I really enjoyed reading it and would like to congratulate the authors on the nice work.
strength_weakness: strength_weakness: * It was a good choice, nonetheless, to explore non-neural models for the prediction of mechanical properties, which is where the authors have the fewest data points (6K). However, given the limited data on coformers (about 7K co-crystals), I also would have expected a simpler baseline like a random forest for estimating likelihood of co-crystal formation. 
* Clarity could be improved in places, since the authors are working with lots of different representations and models throughout their pipeline. Perhaps this could be better visualized in a way that gives a quick overview of the data type, number of data points, and model architecture, for each model included in the pipeline.
* The model may have benefited from a more advanced hyperparameter tuning for the generation of coformers (e.g., Optuna rather than grid search), as the models and their final performance can be quite sensitive to the selected hyperparameters. For instance, %validity >99% should be attainable in all cases for the molecular generative models presented here, to which the models were close but not all quite there.
* Some of the notation is not well-explained. For instance, in the target coformers equation, what is S? In many of the tables and figures, the axes labels or certain notation is implied, which may be mentioned somewhere else in the paper, but it would improve readability if readers could have a reminder of what that property is (and if it is bounded, a percentage, etc). Also, where there are error bars, it should be stated what these are (e.g., in tables, figures). It is a long paper so this would be appreciated, I found myself flipping back and forth a lot.
* To better understand the diversity of the chemical space spanned by generated molecules, it would have been interesting if the authors better quantified that of the training set as well, and if different training/testing strategies could better assess the generalizability to new coformers rather than simply splitting the data. Were not certain classes of coformers more represented than others?
* Given that the authors want others to use their pipeline to discover new coformers, it would have been great to quantify the coverage of chemical space of the model by metrics other than % validity and % novelty. Did the authors consider such metrics (or even dimensionality reduction visualization/techniques) to get a better sense of the chemical space coverage?
questions: questions: * When predicting co-crystal formation, do the authors also have any negative pairs, e.g., molecules that definitely do not co-crystalize? I would assume the bias is towards positive pairs, which perhaps leads to overestimation of the likelihood of co-crystallization for any new pair of molecules, but this was not fully clear to me.
* In the validation case studies (section 5.4), it is unclear if molecules with overlap or significant molecular similarity to nicorandil, rivaroxaban, and paracetamol were excluded from the training set before using the model for these tasks, otherwise there may be data leakage and the conclusions are not so meaningful.
* I see that molecular descriptors (I am guessing 2D?) were used for representing the coformers for mechanical property prediction. Were the fingerprints also explored here, to less effectiveness? Curious since it seems fingerprints were used in other parts of the pipeline, and not sure if descriptors were necessary for good performance here or not.
* Related to the above question, how important is 3D information? I would expect a lot, and wondering if the authors considered including 3D atomic-level descriptors as well in any of their models.
limitations: limitations: * In addition to the final median probability scores attained by each model for the molecular optimization tasks explored in this study, it would be relevant and interesting to report the sample efficiency. This is a much more relevant metric in this context of molecular optimization, rather than the final probability score, as it gives a measure of how quickly a model learns (e.g., how many oracle calls are needed to obtain a specific score). Recommend to include an estimate of the sample efficiency in a revised version. 
* Is the “median probability score” not in fact a likelihood of sampling specific coformers (and not really a probability)? If so (or if not), I think this could be presented more clearly, as what exactly this metric is was not clear to me. It was probably defined somewhere in the paper but I could not find it easily.
* One big limitation is that it would have been relevant to compare the models for the coformer optimization tasks to a simple baseline, such as a “virtual screening” of the coformers. The advantage of using the generative model presented herein should be that better molecules, and thus better scores, should be achievable with the generative model than by simply screening the database (given a fixed sample budget), but this has not been demonstrated in the current study. It would be a pretty simple/standard baseline since the reference database is small (~7K coformers).
paper_summary: paper_summary: This paper presents GEMCODE, the first co-crystal design AI pipeline, which consists of four components:
- SMILES-based models for coformer generation.
- Classification models for co-crystal property prediction.
- An evolutionary algorithm for coformer optimization.
- A GNN for prediction of the probability of co-crystal formation. 

In addition, experiments are carried out on each component and the entire pipeline.
strength_weakness: strength_weakness: 1. The experimental results of property prediction show that the three mechanical properties are inherently elusive, because the accuracy and F1 score of most models in Figure 2 are less than 0.8, which is not ideal for binary classification tasks. Therefore, my concern is that this paper is a pioneer in both the property prediction and generation/optimization of co-crystals, and will the defects (or biases) of the property prediction models be brought into the generation/optimization?

2. In GEMCODE, coformer de novo generation and optimization are divided into two steps, which I think is to follow the tradition in the field of drug discovery. However, in drug discovery, the main purpose of de novo generation is to find molecular candidates whose "main properties" (such as docking scores) are satisfactory, while the main purpose of optimization is to optimize some other properties (such as toxicity and solubility) while maintaining the "main properties" (usually by a similarity constraint). The property objectives of de novo generation and optimization in GEMCODE is the same (without any similarity constraint), which makes me doubt the rationality of the pipeline.

3. Admittedly, I am not familiar with the co-crystal's chemical background, so it is difficult for me to evaluate GEMCODE's chemical validity. Actually, I also did not find any paper on co-crystals in top conferences on machine learning, which makes me question this paper's suitability for publication at this conference. Overall, the main contribution of this paper is to model the computational pipeline of co-crystal design, and to apply existing AI techniques to it. Therefore, in terms of style, I think the chemistry and cheminformatics community may be more appropriate for this paper than a machine learning conference.
questions: questions: 1. In Table 1, the "diversity of targets" are all greater than 0.9. I'm not sure if they are values for internal diversity, and if they are, the diversity values are quite large for sets containing hundreds of molecules. From a molecular design perspective, this suggests that the design objective may be so simple that many very dissimilar molecules fit it.

2. In the paper, the tabletability of co-crystals is presented as a target to be improved. However, in GEMCODE, the three properties related to tabletability are all binary variables, so for a coformer molecule, its tabletability is also represented by one binary variable. Therefore, is the word "enhance" inappropriate for tabletability?

3. Of the three coformer molecules generated in Table 3, two contain ions. I'm not sure if ions are common in coformer molecules, but I think they're not common in the ChEMBL database. So, is there a significant distribution difference between the pre-trained and fine-tuned datasets?
limitations: limitations: Nothing beyond what has already been stated in Section 6.

Comments: title: Reviewer Response
Reviewer Response
References used in the rebuttal
comment: Dear Reviewers and Area Chairs,

As the end of the discussion period is approaching, we would like to say a few words. 

We extend our sincere gratitude to the reviewers. __Thank you very much for your work!__ Your insightful comments and constructive criticism prompted us to significantly improve the quality of our paper. As mentioned before, we will provide additional clarifications, discussions, as well as new experimental results supporting our original claims in the camera-ready submission.

We are very happy to have addressed most (if not all) of the concerns raised in the initial reviews, which resulted in __3/4 reviewers recommending to accept our work__ for publication with an __average rating of 6.0__. We would like to thank the reviewers also for actively engaging in the discussions!

Notably, regardless of the initial review ratings, __all the reviewers acknowledged the novelty and significance of our work__. We are very pleased to have received such evaluations and kindly ask the area chairs and meta-reviewers to take this into account while making the final decision on the paper. A few quotes highlighting this point are:  

  - _"The manuscript tackles a very interesting and impactful problem in pharmaceutical production, and overall a long-standing challenge in crystal engineering. To my knowledge, this is one of the first efforts to use generative ML to design co-formers..."_
  - _"...this is a novel task that the authors have decided to tackle with GEMCODE."_
  - _"This paper highlights an important application of generative AI to a new domain."_  
  - _"This paper is the first to introduce generative AI into co-crystal design, which is a very important topic for drug development and other fields."_
  - _"...this paper is a pioneer in both the property prediction and generation/optimization of co-crystals..."_

We keep working on the refinements to GEMCODE and look forward to the author notification!

Best regards,  
Authors
Thank you for the clarification. I can understand how the moderate size of the co-crystal dataset would make a scaffold split challenging.

Splitting by Tanimoto similarity is a good idea, and will be a good experiment to do in addition to a random split, to better assess OOD generalization.
We greatly appreciate the reviewer’s thoughtful comments on our paper.  
Below, we provide additional considerations and experimental results to resolve the remaining issues:

1. We acknowledge that there were inaccuracies in our response. Undoubtedly, polymorphism of cocrystals is a significant issue. We will include statements emphasizing its importance for pharmaceuticals. However, since our approach is data-driven, and the data on polymorphism is limited, it was impossible to incorporate it into the predictive models in a meaningful way. We will update the Limitations section accordingly. Thanks to the reviewer, we now identify the influence of polymorphism as a prospective research direction.

2. We fully agree with the reviewer that, in general, there are a number of crystal parameters to consider. However, since we are looking for co-crystals in the field of pharmaceuticals, directional interactions other than hydrogen bonds (such as halogen bonds, chalcogen bonds, pnictogen bonds) are rare. Considering the full set of such interactions is beyond the scope of the study. Certainly, for better generalization of the approach, the set of descriptors for the interplanar criterion should be further improved, but this requires additional data (i.e., those co-crystals having halogen, chalcogen bonds, etc.), currently not available. Although GEMCODE is a stable solution for pharmaceutical co-crystals, the pipeline may not fit so well for other co-crystals, which we will discuss in the updated Limitations section. 

3. In response, we plotted distributions of Tanimoto similarity between the molecules generated by the GAN (https://anonymous.4open.science/r/GEMCODE/rebuttal/self_histograms/GAN2.png), and the transformer-based VAE (https://anonymous.4open.science/r/GEMCODE/rebuttal/self_histograms/VAE2.png) and CVAE (https://anonymous.4open.science/r/GEMCODE/rebuttal/self_histograms/CVAE2.png). For each model, the average Tanimoto similarity is between 0.70 and 0.75. On one hand, this highlights sufficient diversity of the molecules. On the other hand, relatively high average similarity was expected, since all generated coformers relate to the same drug. The latter fact aligns well with the observation that the distribution of CVAE is shifted towards 1 due to the “condition” block of the architecture enforcing the target properties of coformers specific to the drug.  
We thank the reviewer for suggesting this analysis. We find it very useful, as it supports our key findings but provides additional insights. We acknowledge a moderate fraction of “duplicates” with Tanimoto similarity = 1 among the predicted molecules. Dropping such molecules brings down the percent of target coformers of Table 1 to 2.23%, 1.68% and 5.63% for GAN, T-VAE and T-CVAE, respectively. Please note that these adjustments do not qualitatively change our results and conclusions. We will update the corresponding sections of the paper accordingly.

4. We agree that baselines from Guacamol (https://github.com/BenevolentAI/guacamol_baselines) can be used to optimize any molecule in SMILES notation with a given objective. However, unlike the Guacamol tasks, the co-crystal design task is multi-objective. The algorithms from Guacamol_baselines (e.g., the noted GraphGA) are focused on single-objective tasks. There are papers where existing Guacamol tasks are used in the mutli-objective formulation, but they require different optimization techniques [1].  
Nevertheless, __we developed a multi-objective modification of GraphGA__ (with Pareto dominance based fitness; the implementation is available at https://anonymous.4open.science/r/GEMCODE/GraphGA_baseline/graphga.py), suitable for co-crystal design tasks. In the __new experiments__, we started from a random subset of co-crystals (with the same population size and number of iterations as was used in GEMOL). GraphGA demonstrated inferior results to GEMOL in terms of target mechanical properties (see plots at https://anonymous.4open.science/r/GEMCODE/GraphGA_baseline/GEMOL_vs_GraphGA.png). More specifically, the highest average probability over all runs was 0.94 (GraphGA) vs 0.95 (GEMOL) for unobstructed planes, 0.63 vs 0.72 for orthogonal planes and 0.12 vs 0.18 for h-bonds bridging. In addition, the convergence of GraphGA was unstable. Ultimately, the hybrid approach of GEMCODE resulted in 21.3% of target molecules on average, while GraphGA produced 20.5% of target molecules of lower quality in terms of predicted mechanical properties. Therefore, we conclude that GEMCODE clearly outperforms the multi-objective version of GraphGA. We will add an appendix section to the camera-ready submission to highlight this comparison. Also, we will further investigate options to extend the set of baselines.

In light of all this, we kindly ask the reviewer to update the initial review.

[1] Optimized drug design using multiobjective evolutionary algorithms with SELFIES //arXiv preprint arXiv:2405.00401. - 2024.
Thank you for your detailed rebuttal and for addressing the concerns raised. I decide to raise my score up to 5.
I very much appreciate the authors for the comprehensive rebuttal and the additional plots. Unfortunately, my key reservations have not been solved:

1. While the authors suggest that the influence of polymorphism on predicting mechanical properties is insignificant due to its low reporting percentage, I must respectfully disagree. Polymorphism is very much underreported because chemists often do not explore further once a crystal structure is identified (especially in non-pharmaceutical contexts). In crystal structure prediction literature, hundreds of local energy minima are almost always predicted for any given compound. Often the case, these predictions are later validated experimentally. 
2. For the correlation between crystal structures' mechanical properties and tableatability, can the authors explain how the diversity of intermolecular interaction affects the prediction? Here, the authors only consider a few criteria, but I am very certain crystals with e.g. halogen bonds would exhibit no hydrogen bonds between planes (hence good score here) but in reality, act very similar to hydrogen bonds. The cited CrystEngComm paper analyzed only 30 crystals and indeed claimed 'While this tool was not intended to be used alone as a predictor of mechanical properties, it _seems to correlate well_ with mechanical properties'. For reasons here and the point above, I do think the claim of tableatability needs to be validated with wet lab results (and hence more suitable for another venue).
3. In Section 5.4, the authors highlight the generation of three new coformers, all of which are carboxylic acids. How is the diversity (in terms of Tanimoto similarity) between generated structures (other than novelty, and I appreciate the histograms)?
4. Maybe this is a misunderstanding, but if the models predict the scores given a pair of SMILES string, why can the authors not employ typical baselines such as Guacamol? Sure, the search space would have to be different and you cannot retrain your scoring function, but GraphGA can surely yield molecules with optimized scores.
We are very grateful to the reviewer for consideration and the decision to increase the score!
In the author console, we still see the rating of 6 though. It could be due to a technical issue with OpenReview. Please make sure to have edited and submitted the updated rating in the original review ("Edit" > "Official Review" > ... > "Submit"), such that it is reflected in the average rating of our work. Thank you!
I thank the authors for their comprehensive rebuttal. The work done here is comprehensive and well motivated. And the author response has also sufficiently addressed my concerns. 

I believe the inclusion of some sort of theoretical validation of any future novel co-crystals would be highly valuable, if at all possible. Experimental validation is expensive, and the authors have demonstrated that previously reported literature results not in the validation set confirm that GEMCODE is effective. When open-sourced, having some sort of way to quickly provide theoretical estimates would be very useful for potential users; almost providing a virtual screening sort of pipeline. Perhaps even a regression model fitted to the available data to act as a proxy.

I am willing to **increase review score to 7**.
We appreciate the reviewer's prompt reply and the chance to provide additional clarifications. We also fully agree that one should not neglect the importance of data splitting and sufficient empirical evidence supporting such study design choices should be provided. We are aware of molecular scaffolding as a method commonly used in drug design to evaluate out-of-distribution generalization capability of predictive models. However, adapting this method to our co-crystal dataset presents several challenges. First, each sample consists of two coformer molecules with different scaffolds (https://anonymous.4open.science/r/GEMCODE/rebuttal/cocrystals.png). Second, the coformers can exhibit a variety of structures, resulting in a large number of selected scaffolds for analysis (approximately 1000 scaffolds based on our preliminary experiments using the Murcko decomposition method). In contrast, in drug design applications (https://greglandrum.github.io/rdkit-blog/posts/2024-05-31-scaffold-splits-and-murcko-scaffolds1.html), researchers typically work with specific compound classes where it is feasible to identify a relatively small number of scaffolds (10-20). Therefore, given the moderate size of the co-crystal dataset, scaffolding might not produce statistically significant evaluations.

At the moment, we identify another  approach as the most promising to further shed light on out-of-distribution generalization. In essence, the data can be split based on Tanimoto Similarity, maximizing dissimilarity between the molecules of different subsets. We are going to conduct such experiments in-depth and include the key findings to the Appendix of the camera-ready submission. We are thankful to the reviewer for rigorous evaluation of our work and attention to detail.
Thank you to the authors for the detailed response and for the clarifications.

One thing that I disagree with still is the comment that
> To our knowledge, there is no common stratification strategy for splitting the data in this domain... Based on our empirical results, we do not expect any significant change in results for alternative data splits.
I do not think this is something that should be swept under the rug. Doing a random split is a standard way of assessing in-distribution generalization, but can you think of other ways to split the data to asses out-of-distribution generalization? I think there are a few interesting experiments you could do, which are not difficult, that would add immensely to the value and utillity of the work. These are also standard experiments to run in applied ML papers.

Related to the above point, there also seems to in general be a lot of confusion about what properties/features are spanned by the data (looking at the other reviewer comments as well as mine), which I think could be better clarified.
__References:__

[1] Heng, T., Yang, D., Wang, R., Zhang, L., Lu, Y., & Du, G. (2021). Progress in research on Artificial Intelligence applied to polymorphism and cocrystal prediction. ACS omega, 6(24), 15543-15550.  

[2] Bryant, M. J., Maloney, A. G. P., & Sykes, R. A. (2018). Predicting mechanical properties of crystalline materials through topological analysis. CrystEngComm, 20(19), 2698-2704.  

[3] Mannava, M. C., Gunnam, A., Lodagekar, A., Shastri, N. R., Nangia, A. K., & Solomon, K. A. (2021). Enhanced solubility, permeability, and tabletability of nicorandil by salt and cocrystal formation. CrystEngComm, 23(1), 227-237.  

[4] Tripp, Austin, and José Miguel Hernández-Lobato. "Genetic algorithms are strong baselines for molecule generation." arXiv preprint arXiv:2310.09267 (2023).  

[5] Ye Z. H. et al. Searching new cocrystal structures of CL-20 and HMX via evolutionary algorithm and machine learning potential //Journal of Materials Informatics. – 2024. – Т. 4. – №. 2.  

[6] Brown N. et al. GuacaMol: benchmarking models for de novo molecular design //Journal of chemical information and modeling. – 2019. – Т. 59. – №. 3. – С. 1096-1108.  

[7] Thomas M. et al. Augmented Hill-Climb increases reinforcement learning efficiency for language-based de novo molecule generation //Journal of cheminformatics. – 2022. – Т. 14. – №. 1. – С. 68.  

[8] Birolo, R., Özçelik, R., Aramini, A., Gobetto, R., Chierotti, M. R., & Grisoni, F. (2024). Deep Supramolecular Language Processing for Co-crystal Prediction. (https://openreview.net/forum?id=bQ9d2hzjW4&noteId=9SoErgR0kb)