Multi-objective de novo molecular design of organic structure-directing agents for zeolites using nature-inspired ant colony optimization
Organic structure-directing agents (OSDAs) are often employed for synthesis of zeolites with desired frameworks. A priori prediction of such OSDAs has mainly relied on the interaction energies between OSDAs and zeolite frameworks, without cost considerations. For practical purposes, the cost of OSDAs becomes a critical issue. Therefore, the development of a computational de novo prediction methodology that can speed up the trial-and-error cycle in search for less expensive OSDAs is desired. This study utilized a nature-inspired ant colony optimization method to predict physicochemically and/or economically preferable OSDAs, while also taking molecular similarity and heuristics of zeolite synthesis into consideration. The prediction results included experimentally known OSDAs, candidates having structures closely related to known OSDAs, and novel ones, suggesting the applicability of this approach.
Introduction<!>Results and Discussion<!>Conclusions<!>Computational Methods<!>Conflicts of interest
<p>The deployment of materials into commercial products needs to satisfy multiple requirements including properties, durability, economic efficiency, and environmental regulations. A class of microporous silicate-based materials known as zeolites has overcome the valley of death to be adapted in several applications such as catalysis and separation. [1] Zeolites can be synthesized under hydrothermal conditions, where inorganic and/or organic cations, together with water molecules, fit within the interior spaces of the anionic silicate-based frameworks to direct the crystallization of specific crystalline phases. [2] Diversification and tuning of organic structure-directing agents (OSDAs), which are typically organic cations, have been a major driving force in expanding the number of available zeolite structures to over 240. [3,4] The use of OSDAs can, however, increase the production cost of zeolites and entails additional processes for removal of OSDAs (mostly combustion) and wastewater treatment. In fact, industrially important zeolites such as FAU, FER, MOR, LTA and LTL zeolites can be produced without using OSDAs. [5] On the other hand, structures including *BEA, CHA, AEI, and CON zeolites, which typically need OSDAs for their syntheses, have also gathered interest due to their remarkable catalytic properties. [6][7][8] In order to advance the implementation of such cost-intensive materials, modification of typical synthetic routes is necessary. A previous research proposed three strategies for the cost reduction of synthesis of zeolites with OSDAs, namely, OSDA-free synthesis, recycle of OSDAs, and replacement of OSDAs. [9] Conventional OSDA-free synthesis has been explored for decades in a wide range of synthesis conditions to produce about 30 aluminosilicate zeolites. [5] The recent expansion of OSDA-free synthetic routes has mainly been driven by the use of seed crystals, which are synthesized from conventional routes using OSDAs. [10] The secondary use of the synthesized product as seeds can substantially reduce the amount of OSDAs required for synthesis, [11] although the number of applicable systems seems to be limited. In the second strategy, OSDAs in zeolites are recycled by breaking and reuniting cyclic ketal groups in OSDA molecules. [12] However, the implementation of recyclable functional groups significantly narrows the freedom of OSDA design, limiting the applicability of the approach. Replacement of complex and expensive OSDAs with simpler ones is successful in several cases. [9,[13][14][15][16][17] Furthermore, instead of relying on time-consuming trialand-error experiments to find reasonable alternative OSDAs, molecular modeling can tackle the problem. OSDA's ability to direct the crystallization of zeolites is well understood to be mainly derived from van der Waals interactions and charge compensation between OSDAs and inorganic zeolite frameworks. Even though crystallization of zeolites is strongly governed by kinetics, [18] computational modeling of the host-guest interactions can accurately describe the outcomes from hydrothermal synthesis, [19][20][21][22] selectivity against competing phases, [23,24] configurations of OSDAs, [25][26][27] and location of Al in zeolites with OSDAs. [28] Employment of the host-guest interactions as the objective function has been realized in de novo molecular design methodologies for zeolite-OSDA systems. The first example of the methodology demonstrated that stochastic algorithms can grow organic molecules inside the zeolite frameworks. [19] However, the growth algorithm also yields organic molecules that are difficult or impossible to be synthesized. This issue is solved in recent studies by Deem and his co-workers, [29] where an evolutionary algorithm treats an array of readily available reagents and known chemical reactions as a genetic sequence for producing candidate OSDAs. The method can predict OSDAs along with synthetic pathways and stabilization energies, which are useful to verify in zeolite synthesis. [16,22,24,30,31] However, this evolutionary algorithm may sometimes predict trivial, unrealistic, and/or economically less advantageous candidates. [29,31] Design of OSDAs is a multi-objective problem. Several issues, including stabilization energy, cost of reagents, synthetic pathways, and environmental compatibility, must be considered simultaneously, making it a complicated Pareto optimization problem. The difficulty is apparent in the network of organic chemistry where molecules are regarded as nodes and chemical reactions are regarded as edges. [32,33] The exploration of the large chemical network to find the optimal compound is computationally demanding and requires meta-heuristic approaches that are particularly suited for network-based problems. The present study solves the multi-objective problem of OSDA design using a nature-inspired ant colony optimization (ACO) algorithm. [34] The methodology can optimize stabilization energy of zeolite-OSDA systems by considering the cost of reagents and molecular similarity, thus screening the large database of available reagents. The stabilization energy and/or cost of the predicted results outperform current state-of-the-art OSDAs, proving the efficiency of the method used.</p><!><p>Although synthesis of OSDAs occasionally requires multiple synthetic steps using various organic reactions, most of the practically achievable OSDAs are produced by alkylation of amines with organic halides. If we focus on particular organic reactions, de novo design algorithms including ACO can be efficiently used. Therefore, the scope of OSDAs considered here is limited to the products of alkylation of commercially available amines and halides. ACO is a general meta-heuristic optimization algorithm inspired by nature, where artificial ants move stochastically within a network while tracking pheromones deposited at each edge. [34] High-pheromone edges are selected with high probability because this means the edges are frequently visited by well-performing ants. The amount of pheromone is determined from scores of an objective function to be optimized. The ability of ACO to find the best route from a network is suited for solving network-based problems such as the traveling salesman problem. [35] ACO has also been applied to chemical networks, for instance, to reductive amination for construction of predictive models of structure-activity relationship. [36,37] The efficiency of ACO in chemical networks was demonstrated in the implementation of the de novo design in automated synthesis and assay systems. Our computational framework involves data collection for the library of organic compounds, [38][39][40][41] synthesis of organic cations using virtual organic reactions, addition of OSDAs into cavities of zeolites, exploration of different configurations of OSDAs occluded in zeolites, computational estimation of stabilization energies by molecular dynamics (MD), identification of the winner ants, and updating the pheromones. Fig. 1 schematically represents the workflow of our computational protocol. Amines are regarded as the first nodes while halides are the terminal nodes to be visited by artificial ants. Each artificial ant collects an amine and a halide through a stochastic process by following high-pheromone nodes. Virtual alkylation reaction produces a quaternary ammonium cation, which is subsequently evaluated using multiple filters. The first filter eliminates cyclic ammoniums with significant distortions, which can be detected by counting the numbers of π orbitals in an ammonium ring. The second filter is a degree of hydrophobicity of cations. It is generally accepted that OSDAs need to have a moderate hydrophobicity to be dissolved in aqueous solutions and to interact with hydrophobic zeolite frameworks. [42] The hydrophobicity of OSDAs can be determined by a partitioning experiment that evaluates the degree of phase transfer of OSDA salts from water to organic solvent, mostly chloroform or octanol. The chloroform-water partition coefficients have been Fig. 1. Schematic illustration of the computational workflow. A colony contains 10 artificial ants where each agent stochastically selects an amine and a halide with the probability proportional to the number of pheromones. The distances between nodes are defined by the cost parameter η and the adjustable parameter β. Virtual alkylation reaction produces quaternary ammonium cation, which is subsequently screened with computationally efficient criteria. The candidate cation that satisfies the criteria is subjected to the MD run to calculate the stabilization energy, which is used to update the number of pheromones for subsequent ants. correlated with C/N + atomic ratios of efficient OSDAs. [43] Following this knowledge, we reject candidates with too low (≤ 5) and too high (≥ 18) C/N + . The rotational freedom of OSDAs is another effective metric to screen candidate molecules. Molecules with many rotatable bonds should be avoided because their conformations can be changed under high temperature in hydrothermal synthesis, which can reduce the selectivity of zeolite crystallization. [43] The number of rotatable bonds is restricted to less than 3 in accordance to the previous report. [43] After removing molecules that are too large to stay in the inner space of a given zeolite, multi-step MD runs are conducted to evaluate the stabilization energy, which becomes the score of the ant. After running 10 ants, the pair of nodes visited by the best ant receives pheromones.  [a] (2.0 [b] ) Tetraethylammonium 14 [a] (0.50 [c] ) Tetrapropylammonium 32 [a] (4.8 [d] )</p><p>[a] Calculated from the cost of amines and halides. Considering the high cost parameter η of 1, we explored the chemical network using ACO to find alternative OSDAs that can direct CHA crystallization. While some of the metaheuristic algorithms are virtually black box, the adaptive process of ACO is easily interpretable by tracking the amount of pheromone deposited on reagents. Fig. 2 shows the trajectory of normalized pheromone amount in the typical first 100 steps for CHA. As can be seen, all reagents started from an equivalent amount of pheromone, which exponentially decayed as ACO optimization progressed.</p><p>After 10 steps, one ant selecting A1 produced a candidate OSDA that outperformed the others selected by other ants. At this point, the candidate OSDA synthesized by A1 was a global optimum that could receive an extra amount of pheromone, fostering the exploration of the other candidates that use A1 as a building block. One of the common straightforward approaches for molecular design is the use of structurally similar molecules.</p><p>To mimic this design strategy, we implemented a function to consider molecular similarity. [54] Amines that are structurally similar to A1, such as A2, A3 and A4, got different amount of extra pheromones after A1 became the global best. With the increased probability as the step proceeded, candidate OSDAs derived from A2, A3 and A4 were selected by ants and won within colonies as shown in the sharp uptakes of the normalized amounts of pheromones (i.e., probabilities) until step 80. At step 80, A5, the amine that had increased its amount of pheromone at step 58, became the global best by reacting with H1, resulting in the decrease of A1-A4 and the increase of A5 and its structurally related amines until step 90, where the combination of A6 and H1 outperformed the previous global best. In the corresponding trajectory for halides, H1 remained the global best, although other halides (H2, H3, and H4) were frequently selected and won local ant colonies. Due to the max-min functionality, repeated winnings of the global best did not show uptakes, although the additional pheromone derived from molecular similarity increased the probability to be selected. Thus, ACO can sample candidate OSDAs both similar and different from known ones. [37] Fig. 3. OSDAs for CHA computationally predicted without considering the cost of chemicals (a cost adjustable parameter, β = 0) and their corresponding stabilization energies Es and cost parameters η. An asterisk denotes the experimentally proven OSDAs for CHA. Fig. 3 summarizes some of the computationally predicted OSDAs for CHA. The automated filling algorithm introduced three molecules per unit cell in most cases, which is consistent with experimental observations [53] and the number of cages in CHA. [3] Obviously, the ACO algorithm can propose candidate OSDAs comparable with the experimentally proven OSDAs, summarized in Table S1, in terms of stabilization energies. The fact that ACO identified experimentally proven OSDAs for CHA such as 1 suggests the excellent predicting ability of the computational protocol. The most stabilized structure of the zeolite-OSDA complex shown in Fig. 4a corroborates the strong stabilization of the system where symmetric 1 structurally fits with the cylindrical cage of CHA. OSDA 4 is structurally related to benzyltrimethylammonium, which can be used to reduce the amount of 1 required for the synthesis of CHA. [9] 2 and 5 are slightly smaller than an experimentally proven OSDA, 34 (see Table S1). The most stabilized complex of CHA and 2 shown in Fig. 4b indicates that the candidate OSDA is tightly confined in the cage of CHA. The only difference between predicted 8 and known 37 (see Table S1) is the two bridging carbons, which makes 37 slightly larger. As was implemented, but importantly, all predicted candidates had moderate hydrophobicity and rigid frameworks in terms of C/N + and the number of rotatable bonds, which is essential for successful crystallization of zeolites as explained earlier. [43] The stabilization energy of 3 was comparable with experimentally known OSDAs and its cost parameter was significantly more advantageous than 1, the typically used OSDA for CHA, suggesting the possibility to improve economic efficiency. The current methodology is applicable to any zeolites as long as the structure has an inner cavity to host OSDAs. To demonstrate this, we applied the de novo OSDA design framework to another small-pore zeolite, AEI, which is structurally related to CHA. [3] Compared to the cylindrical nature of the cage of CHA, the cage for AEI is an unsymmetrical avocado-like structure as shown in Fig. 4c  and d, which is one of the reasons why different OSDAs are required for the crystallization of CHA and AEI. As summarized in Table S2, several OSDAs have been used to direct the crystallization of AEI by filling the distinct cage shown in Fig. 4c and d. According to a previous study, the stabilization energies of experimentally proven OSDAs for AEI ranges from −13.2 to −16.9 kJ/molSi. [31] The ACO algorithm found comparable candidate OSDAs as summarized in Fig. 5. Among them, 9 outperformed the stabilization energy of state-of-the-art experimentally proven OSDAs. [31] The unsymmetrical molecular structure of 9 seems to be appropriate for the avocado-like cage of AEI as shown in Fig. 4d. ACO also found experimentally proven OSDAs for AEI such as 16, supporting the broad applicability of the present computational protocol. An illustrative outcome of the functionality to consider molecular similarity was 13, 14, and 15 in that they are similar to experimentally known OSDAs such as 16, 50, and 51 (see Table S2). A notable issue of candidate OSDAs predicted in Fig. 3 and 5 is the high cost parameter η, which can be problematic once the zeolites are subject to large-scale experiments or industrial applications. In the case of CHA, even though some of the candidates including 3 and 4 are economically more efficient than 1, other candidates such as 2, 6, and 7 show prohibitively high cost parameters (see Fig. 3). This can cause convergence to physicochemically appropriate, but economically unrealistic candidates. In order to consider the cost of chemicals, we return to the fundamental equation to run ACO. Previous studies employing ACO for adaptive molecular design used identical "distance" between all reagents. [36,37] As depicted in Fig. 1, our current study uses the cost parameter η as the "distance" to increase the probability that an ant selects lower cost reagents upon optimizing the stabilization energies. This distance can be adjusted using another parameter, β. When β is set to 0, ACO runs without cost consideration (i.e., uniform distance), as shown above in Fig. 3 and 5. Conversely, the cost parameter fully represents the "distance" when β = 1 (see details in Computational Methods). Fig. 6 shows the mapping of candidate OSDAs predicted by a typical run of ACO for AEI with different β, to control the degree of consideration of the cost parameter. With β = 0, artificial ants sampled both economically reasonable and expensive chemicals to find appropriate candidates based solely on the host-guest interactions. As a result, the prediction produced candidates that can be prohibitively too expensive to become alternatives (vide supra). A slight consideration of the cost parameter with larger β reduced the frequency of expensive reagents, as shown in Fig. 6b. This results in less negative stabilization energies, likely due to less chances to sample expensive but physicochemically preferable chemicals. Further increase of β highly limited the chance to sample cost intensive reagents as evident from the scarcity of samplings with η > 500 (Fig. 6c and  6d).</p><p>A previous study pointed out that the restriction of applicable reactions sometimes yields superior molecules compared to predictions with more freedom. [29] This seemed to hold for a case of β = 0.5, where 17 marked stabilization energy of −18.7 kJ/molSi (Fig. 7), outperforming all candidate OSDAs predicted with lower β (Fig. 5) and all experimentally proven OSDAs for AEI (Table S2). Further increase of β to 1 limited the possibility to select the cost intensive reagents more tightly, which resulted in less negative stabilization energy, except for 18 shown in Fig. 7 which defeats all experimentally proven</p><p>OSDAs for AEI (Table S2) in terms of the host-guest interaction. The cost parameter η of the predicted OSDAs with β = 1 was notably smaller than the candidates identified without considering the economic efficiency, while the stabilization energies were comparable to experimentally proven OSDAs and predicted OSDAs without cost consideration. Notably, 19 predicted with β = 1 is one of the experimentally proven OSDAs for AEI that is physicochemically and economically preferable (Fig. 7). Accordingly, our proposed algorithm traced experimental efforts to explore physicochemically and economically reasonable OSDAs for the synthesis of zeolites by tuning the cost parameter β. At cost parameter β = 0.5, ACO generated some experimentally proven OSDAs for CHA, including 20 [14] and 21 [10] (Fig. 7). However, 20 and 21 require another OSDA such as 1 for successful crystallization of CHA, presumably owing to their weaker structure-directing abilities [10,14] as indicated by their less negative stabilization energies. Despite this, 20 and 21 are significantly more advantageous than 1 in terms of economic efficiency as is evident from their lower η values. The function to consider molecular similarity was effective again for the current case, where the algorithm predicted structurally similar but more physicochemically preferable candidate OSDAs including 22, 23, and 24 in the same run.</p><p>Another zeolite that suffers from the high cost of OSDAs is CON, [17] which is a three-dimensional medium-and largepore zeolite with excellent catalytic performance in methanol-to-olefin reactions. [8] The OSDA that can crystallize pure CON zeolite, 25 [61] has a cost parameter of 310 (Fig. 8), which is extremely expensive. Unlike the systems where isolated OSDAs are geometrically fitted within the cages as shown in Fig. 4, more than one OSDA seems to fill the intersection of the multidimensional channels in CON. [17,56] We ran systematic MD simulations for 25, finding that four molecules can be accommodated in the unit cell of CON (see Table S3), which is consistent with the previous study. [56] Fig. 9a shows that a pair of 25 interacted with each other by van der Waals interactions. The ACO run for CON with β = 0 successfully predicted 25 (Fig. 8). With larger β of 0.5, more reasonable candidates were sampled. In the cases of 26, 27, and 28, four molecules were filled in the unit cell of CON to stabilize the zeolite-OSDA complexes. It is known that π-π interactions between OSDAs can create supramolecular self-assemblies to fill large cavities in porous materials. [57] The π systems in 27 and 28 have a potentially positive impact on creating molecular aggregates in inorganic cavities to direct the crystallization of CON. Smaller molecules such as 29 and 30 also showed superior stabilization energies compared to the experimentally proven 25. Six molecules filled the unit cell of CON in the cases of 29 and 30 (Fig. 9b). Very recently, our group reported the successful crystallization of CON in the presence of tetraethylammonium and seed crystals. [17] The only difference between the predicted OSDAs (29 and 30) and tetraethylammonium is the C-C bonds that connect ethyl groups at the terminals. Since the current methodology rejects molecules like tetraethylammonium owing to the high rotational freedom, very similar but more rigid molecules such as 29 and 30 were predicted as the alternative candidate OSDAs for CON, which are promising from both physicochemical and economic viewpoints. Although the current study assumed that sufficiently negative stabilization energy is the indicator of promising structure direction properties, competitions with other zeolite phases can be remarkably critical. [16,23,24,29] In such situations, the mixing of different OSDAs (i.e., the use of more than one OSDA) [10] and the use of seed crystals [17] can lead to successful crystallization, which may economically outperform the existing preparation protocols. Another metric to evaluate the cost of organic synthesis is the number of steps to achieve the desired molecule. [16] Although the current methodology focused only on the alkylation reaction, other reactions with multiple steps can feasibly be considered by using a larger network of organic chemistry and potential "distance" information stored on each node (chemical) and each edge (reaction), including cost, safety and yield, [33] which offers the chance to optimize the production of zeolites from multiple aspects.</p><!><p>Organic structure-directing agents (OSDAs) have massively contributed to the development of zeolite science as they seem to accelerate the discovery of new zeolite materials and to optimize existing zeolites for certain applications. However, these organic molecules can become an economic burden once zeolites are subject to large-scale experiments and/or industrial applications as quantified by the cost parameters. Although computational methods can relieve the cost of screening a large number of possible molecules, it can suggest unfeasible, prohibitively expensive, unrealistic, and/or trivial candidates. To design OSDAs, this study employed a nature-inspired ant colony optimization (ACO) de novo design algorithm with implementation of chemical heuristics and molecular similarity to construct the design framework for optimization of physicochemical objective function with consideration of economic efficiency. The predicted candidates contain experimentally proven OSDAs, molecules similar to experimentally proven ones, and novel molecules owing to the unique advantage of ACO molecular design of simultaneously sampling structural analogues of known molecules and completely new ones. [37] Several recent studies have revealed that OSDAs can control morphology, [25] chirality, [24] and atomic configurations [28] of zeolites, which can be evaluated by computational methodologies. It is therefore anticipated that multi-objective design of OSDAs to simultaneously optimize these zeolite properties can be approached by ACO in near future.</p><!><p>The library of available organic compounds was retrieved from commercial suppliers. [38][39][40][41] This dataset containing 97,132 compounds was reduced by eliminating molecules with elements other than carbon, hydrogen, nitrogen, chlorine, bromine, and iodine. Molecules unable to participate in the alkylation reaction were taken out of consideration. The prices and the purity of chemical compounds were retrieved from Aldrich Market Select and websites accessed from Japan. [38] The cost parameter, η, was estimated from an empirical model described in a previous study [32] and multiplied by the purity and molar weight of the corresponding chemicals. Virtual reactions were encoded using MarvinSketch [58] and conducted until saturation. The resulting molecule was structurally optimized using a universal force field. [59] The volume, [60] the number of π orbitals, and the number of rotatable bonds of a molecule were calculated using RDKit software. [61] Crystallographic information of idealized zeolites was retrieved from a database provided by the International Zeolite Association. [3] The locations of molecules in zeolite pores were identified using our in-house computational protocol. First, the Voronoi network in zeolite pores was constructed by Zeo++ software. [62] Among nodes of the Voronoi network, the location with least steric hindrance for the molecule was identified by counting the number of atoms too close to each other. After finding the best location, the molecule was randomly rotated for 100 different angles to minimize the steric hindrance. Following the previous study, [29] structure optimization and multi-step molecular dynamics (MD) simulation were performed on GULP [63] under periodic boundary condition. Charge-less Dreiding force field [64] was applied to consider the rapid motion of organic molecules [2] and vibration of zeolite frameworks with 343 K on an NVT ensemble. If the connectivity of the model was changed during the MD run, the candidate OSDA was rejected. An empty zeolite and a free molecule were also subjected to the MD run with the same setting except that the molecule was simulated as a discrete model. The values of the last 5 ps of the final MD run with the production time of 30 ps were averaged to calculate representative energies. [29] Stabilization energy 𝐸 " was computed from the following equation:</p><p>𝐸 " = 𝐸 $%&'()* − 𝑛𝐸 -./0 − 𝐸 1)% (23)   where 𝐸 $%&'()* is the energy of the zeolite-OSDA complex, 𝑛 is the number of OSDAs, 𝐸 -./0 is the energy of the free OSDA, and 𝐸 1)%( 23) is the energy of the empty zeolite. The composition of zeolites was assumed to be purely siliceous throughout the computations. Energies were normalized by the number of silicon in zeolites. Some of the values were used from a previous study. [31] In ACO, an artificial ant at node 𝑖 selects the next node 𝑗 with the following probability: where 𝜏 78 is the amount of pheromone deposited on the between 𝑖 and 𝑗, 𝜂 is the closeness between 𝑖 and 𝑗, and 𝛽 is an adjustable parameter. The scoring function 𝛥𝜏 78 derived from the 𝑖 th amine and the 𝑗 th halide to be maximized was defined as, where 𝑡 is the ACO time step and ρ = 0.95. The size of ant colony was set to 10. Morgan fingerprint [65] was used with a radius of 4 to calculate a bit vector with 2048 bits. The amount of additional pheromone based on structural similarity between molecular structure was computed by the following equation:</p><p>where 𝑠 is the Tanimoto similarity index [54] between a reagent and the global best at that moment. Max-min filter [66] was applied for all pheromones. The maximum amount of pheromone 𝜏 &D* and the minimum amount of pheromone 𝜏 &2F were calculated according to the following equations: where 𝜏 `(%ED( E)"3 was the amount of pheromone of the global best and 𝑛 was the number of amines or the number of halides. The amount of pheromone was updated when at least one ant among a colony could calculate the stabilization energy from MD.</p><!><p>There are no conflicts to declare.</p>