Computer-aided molecular design: An introduction and review of tools, applications, and solution techniques

TL;DR

This paper reviews computer-aided molecular design (CAMD), covering its methods, optimization techniques, and applications, aiming to serve both beginners and experts in the field.

Contribution

It provides a comprehensive overview of CAMD tools, methods, and applications, including recent optimization strategies and structural feasibility constraints.

Findings

Effective use of QSPRs in CAMD problems

Various optimization techniques for CAMD

Diverse applications of CAMD in molecular design

Abstract

This article provides an introduction to and review of the field of computer-aided molecular design (CAMD). It is intended to be approachable for the absolute beginner as well as useful to the seasoned CAMD practitioner. We begin by discussing various quantitative structure-property relationships (QSPRs) which have been demonstrated to work well with CAMD problems. The methods discussed in this article are (1) group contribution methods, (2) topological indices, and (3) signature descriptors. Next, we present general optimization formulations for various forms of the CAMD problem. Common design constraints are discussed and structural feasibility constraints are provided for the three types of QSPRs addressed. We then detail useful techniques for approaching CAMD optimization problems, including decomposition methods, heuristic approaches, and mathematical programming strategies. Finally, we discuss many applications that have been addressed using CAMD.

Tables

Table 1. Table 1: Sample of available GC methods for predicting various properties of pure compounds

Property	GC methods
Aqueous solubility	Marrero et al. [114], Klopman [100]
Boiling point	Joback and Reid [78], Stein et al. [157], Nannoolal et al. [122], Marrero et al. [113]
Bond dissociation energy	Benson et al. [10]
Critical pressure	Jalowka et al. [76], Joback and Reid [78], Klincewicz et al. [99], Nannoolal et al. [123], Marrero et al. [113]
Critical temperature	Jalowka et al. [76], Joback and Reid [78], Klincewicz et al. [99], Nannoolal et al. [123], Marrero et al. [113]
Critical volume	Klincewicz et al. [99], Nannoolal et al. [123], Marrero et al. [113]
Enthalpy of formation	Cohen et al. [31], Benson [9], Domalski et al. [36], Roganov et al. [144], Joback and Reid [78], Marrero et al. [113]
Enthalpy of fusion	Joback and Reid [78], Marrero et al. [113]
Enthalpy of vaporization	Roganov et al. [144], Joback and Reid [78], Marrero et al. [113], Ceriani et al. [24]
${\text{LC}}_{\text{50}}$	Martin et al. [115] (fathead minnow)
Melting point	Joback and Reid [78], Marrero et al. [113]
Gibbs energy of formation	Cohen et al. [31], Benson [9], Domalski et al. [36], Roganov et al. [144], Joback and Reid [78], Marrero et al. [113]
Heat capacity	Benson et al. [11], Benson [9], Domalski et al. [36], Joback and Reid [78], Kolská et al. [103], Ceriani et al. [24]
Octanol/water partition coefficient	Marrero et al. [114], Platts et al. [140], Klopman et al. [101]
Vapor pressure	Nannoolal et al. [124]
Viscosity	Joback and Reid [78], Sastri et al. [149], Ceriani et al. [23], Cao et al. [22], Nannoolal et al. [125]

Table 2. Table 2: Sample of available TI-based methods for predicting various properties of pure compounds

Property	TI method
Anti-inflammatory activity	Gupta et al. [61], Bajaj [5]
Aqueous diffusion coefficient	Schramke [151]
Aqueous solubility	Hall et al. [65], Katritzky et al. [85]
Biodegradability	Boethling [12]
Blood-brain barrier partition coefficients	Rose et al. [145]
Boiling point	Hall et al. [65], Hosoya [73], Hall et al. [64], Galvez et al. [51]
Critical temperature	Hall et al. [64]
Density	Kier et al. [90], Estrada [43], Katritzky et al. [84]
Enthalpy of formation	Mercader [118]
Enthalpy of fusion	Gharagheizi et al. [59]
Enthalpy of vaporization	Galvez et al. [51]
Flash point	Patel et al. [135]
Heat capacity	Yao et al. [169]
${\text{LC}}_{\text{50}}$(fathead minnow)	Hall et al. [66, 67], Basak et al. [8]
Melting point	Katritzky et al. [84]
Nonspecific local anesthetic activity	Kier et al. [89], Katritzky et al. [84]
Octanol/water partition coefficient	Murray et al. [121]
$\pi$-electron energy (C-C bonds)	Hosoya et al. [74]
Refractive index	Katritzky et al. [84]
Vapor pressure	Katritzky et al. [85]
Viscosity	Kauffman et al. [86]
Water-air partition coefficient	Katritzky et al. [85]

Table 3. Table 3: Summary of CAMD applications and the methodologies used in each case

Application	References
Antigen inhibition activity	Churchwell et al. sd,d,o [29]
Biodiesel additives	Hada et al. gc,d,gt [62]
${\text{CO}}_{\text{2}}$ capture	Gani et al. gc,d,gt,m [55], Bardow et al. sel,d,p [6], Pereira et al. gc,o,p [138], Stavrou et al. sel,d,o,p [156], Burger et al. gc,o,p [18], Lampe et al. gc,d,o,p [105], Gopinath et al. gc,d,o,p [60]
Crystallization solvents	Karunanithi et al. gc,d,gt,p,m [83], Samudra and Sahinidis gc,d,o [147], Austin et al.gc,d,o,p,m [4]
Extractive distillation	Harper et al. gc,d,gt [71], Gani et al. gc,d,gt [55], Papadopoulos and Linke gc,d,h,p [131], van Dyk and Nieuwoudt gc,d,h [40]
Extractive fermentation	Wang and Achenie gc,d,o [164], Papadopoulos and Linke gc,d,h,p [129]
Gas absorption	Odele and Macchietto gc,d,o [127], Buxton et al. gc,d,o,m,p [19], Papadopoulos and Linke gc,d,h,p [130, 131], Bommareddy et al. gc,d,gt,p [14], Zhou et al. gc,d,h,p [173]
HIV-1 protease inhibition activity	Visco et al. sd,d,o [163]
Ionic liquids design	Matsuda et al. gc,gt [116], McLeese et al. ti,h [117], Karunanithi and Mehrkesh gc,d,h [81]
Liquid-liquid extraction	Gani and Brignole gc,d,gt [53], Brignole et al. gc,d,gt [16], Odele and Macchietto gc,d,o [127], Marcoulaki and Kokossis gc,d,h [112], Harper et al. gc,d,gt [71], Harper and Gani gc,d,gt [70], Gani et al. gc,d,gt [55], Karunanithi et al. gc,d,gt,p [82], Austin et al.gc,d,o,p,m [3], Ourique and Telles gc,h [128], Kim and Diwekar gc,d,h,p [93], Papadopoulos and Linke gc,d,h,p [130], Xu and Diwekar gc,d,h [168], Scheffczyk et al. d,o [150], Gebreslassie and Diwekar gc,h [58]
Organic Rankine cycle fluids	Papadopoulos et al. gc,d,h,p [132], Lampe et al. sel,d,p [104], Lampe et al. gc,d,o,p [105]
Pharmaceutical products	Siddhaye et al. ti,o [153]
Polymer design	Venkatasubramanian et al. gc,h [161, 162], Maranas gc,o [111], Vaidyanathan and El-Halwagi gc,o,m [160], Camarda and Maranas ti,o [20], Brown et al. sd,o [17], Eslick et al. ti,h [42], Pavurala and Achenie gc,d,o [136], Zhang et al. gc,o [170]
Reactions solvents	Wang and Achenie gc,d,o [164], Gani et al. sel,d [57], Folić et al. gc,o [48, 49], Struebing et al. gc,d,o [158], Zhou et al. gc,o [171], Austin et al.gc,d,o,m [3], Zhou et al. gc,d,h [172]
Refrigerant design	Joback gc,d,gt [77], Gani et al. gc,d,gt [55], Churi and Achenie gc,d,o [30], Duvedi and Achenie gc,d,o [38, 39], Marcoulaki and Kokossis gc,d,h [112], Sahinidis et al. gc,o [146], Ourique and Telles gc,h [128], Samudra and Sahinidis gc,d,o [148]
Separations (general)	Hostrup et al. gc,sel,d,o,gt,p [75]
Solvents for consumer products and industry	Pistikopoulos and Stefanis gc,o [139], Buxton et al. gc,d,o,m,p [19], Conte et al. gc,d,gt,m [33], Sinha et al. gc,o [154], Sinha et al. sel,o [155], Weis and Visco sd,o [165]
Soybean oil products	Camarda and Sunderesan ti,o [21]
Structural modifications to a fungicide	Raman and Maranas ti,o [141], Chemmangattuvalappil et al. sd,o [28]
Transition metal catalyst design	Chavali et al. ti,h [26], Lin et al. ti,h [106]
VOC recovery	Eden et al. gc,d,gt,p [41]

Equations

$$P=\sum _g{c}_g{n}_g$$

$$P=f\left(\sum _g{c}_g{n}_g\right)$$

$$P=f\left(\sum _{g\in F}{c}_g{n}_g+\sum _{g\in S}{c}_g{n}_g+\sum _{g\in T}{c}_g{n}_g\right)$$

$$I_{g,g^{\prime }}=f_I\left(n_g,n_{g^{\prime }}\right)$$

$$W(G)=1/2\sum _{v,v^{\prime }}d(v,v^{\prime })$$

$$C{I}_E(v,v^{\prime })=\frac{1}{\sqrt{{\delta }_v{\delta }_{v^{\prime }}}}$$

$${}^1\chi =\sum _{\{v,v^{\prime }\}\in E}C{I}_E(v,v^{\prime })=\sum _{\{v,v^{\prime }\}\in E}\frac{1}{\sqrt{{\delta }_v{\delta }_{v^{\prime }}}}$$

$${}^0\chi =\sum _{v\in V}\frac{1}{\sqrt{{\delta }_v}}$$

$${}^i\chi =\sum _{\{v_1,v_2,\dots \}\in V_C^i}\frac{1}{\sqrt{\prod _v{\delta }_v}}$$

$${\delta }_v^{\text{V}}=Z_v^{\text{V}}-h_v$$

$${\delta }_v^{\text{V}}=\frac{Z_v^{\text{V}}-h_v}{Z_v-Z_v^{\text{V}}-1}$$

$${}^1\kappa =\frac{2{(}^1{P}_{\text{max}}){(}^1{P}_{\text{min}})}{{(}^{1}P{)}^2}$$

$${}^2\kappa =\frac{2{(}^2{P}_{\text{max}}){(}^2{P}_{\text{min}})}{{(}^{2}P{)}^2}$$

$${}^3\kappa =\frac{4{(}^3{P}_{\text{max}}){(}^3{P}_{\text{min}})}{{(}^{3}P{)}^2}$$

$${}^0\alpha (v)=c_G(v)$$

$${}^i\alpha (v)=G_S(V^i,E^i,C^i,c_G)$$

$$P=\sum _i\sum _{d\in D^i}{c}_d{}^i{\alpha }_G(d)$$

$$\underset{n}{\min }$$

$$s.t.$$

$$h_1(p,n)\le 0$$

$$h_2(p,n)=0$$

$$s_1(n)\le 0$$

$$s_2(n)=0$$

$$p_k^L\le p_k\le p_k^U\forall k$$

$$n_d^L\le n_d\le n_d^U\forall d$$

$$C(p)=\sum _k{w}_k(p_k^T-p_k{)}^2$$

$$\underset{n,x}{\min }C(n,p,q)$$

$$q=g(x,n,p)$$

$$\sum _i{x}_i=1$$

$$h_1(p,q,n)\le 0$$

$$h_2(p,q,n)=0$$

$$q_j^L\le q_j\le q_j^U$$

$$\underset{n,\mu }{\min }C(n,p,\mu )$$

$$h_1(p,\mu ,n)\le 0$$

$$h_2(p,\mu ,n)=0$$

$${\mu }_w^L\le {\mu }_w\le {\mu }_w^U$$

$$n_d\le n_d^U$$

$$\sum _d{n}_d\le N^U$$

Full text

Keywords: Computer-aided mixture design; Computer-aided molecular design; Integrated product and process design; Group contribution; Topological indices; Signature descriptors

Introduction

The application of chemistry to manipulate the natural world has its earliest examples in metalworking and pottery [134], with some pottery artifacts reported to be as old as 20,000 years [167]. Chemical products have since played an important role in history and left an indelible mark on the way we live and work. Early history included basic incendiary fuels, perfumes, and soap as some of the first widespread uses of chemicals. The modern age has witnessed an unprecedented expansion of chemical products, including pesticides, fuels for transportation and electricity, pharmaceuticals, plastics, and a broad array of industrial and consumer products. True to these historical trends, few things are as pervasive as chemical products in 21st century life, and there are an ever-increasing number of new chemical applications which require specialized compounds.

The process of determining new and suitable chemicals for a certain application can be generally termed chemical product design [34]. Chemical product design has long been a laborious, trial-and-error procedure, limited often by a fixed amount of chemical, time, and financial resources. Design efforts are often high-throughput and tend to focus on a small class of compounds or structural analogues of known chemicals. Accordingly, the so-called “design space”—the set of unique molecular structures considered—is often quite small for these product design problems, especially considering the massive design space of all possible chemical structures. It is then clear that to keep pace with the growing demand for new chemical products and to adequately explore the full chemical design space, other approaches must be considered. Fortunately, the availability and efficiency of computational resources makes these design problems more tenable than ever before. Noteworthy among computational approaches is the field of computer-aided molecular design (CAMD), which leverages the simplicity of semi-empirical quantitative structure-property relationships (QSPRs) in conjunction with fast and efficient numerical optimization algorithms.

CAMD has its roots in the 1980’s, although the general use of computers in chemistry pre-dates this by a few decades. Stated formally, the CAMD problem concerns designing an optimal molecular structure(s) for a certain application. CAMD combines molecular modeling techniques, thermodynamics, and numerical optimization to design good or optimal molecular structures, many of them often completely novel. Advances in chemical modeling in the last few decades have greatly benefited CAMD, and practitioners are now capable of relating chemical structures to properties at several levels of accuracy (molecular mechanics, semi-empirical, ab initio). Though CAMD often uses semi-empirical modeling techniques for their simplicity and efficiency, new approaches incorporating more accurate methods are emerging. Modern combinatorial optimization techniques are also essential for CAMD, enabling the optimization over staggeringly large design spaces which would otherwise be inaccessible (using enumeration algorithms, for example).

This article is intended to provide a description of popular QSPRs used in CAMD, the CAMD problem itself, and several solution approaches to its various forms. There exist several CAMD reviews [52, 126] that provide a thorough coverage of many relevant applications and milestones of CAMD and are good resources for the interested reader. While also providing a review of CAMD literature, the current article first motivates the field from the ground up. In this way, this review is intended to be accessible to readers of all levels as it presents many essential CAMD tools and techniques in detail. We also provide a discussion of various techniques and a review of CAMD applications, so this document will also be of interest to the more advanced CAMD practitioner. We begin in Section 2 by detailing three popular classes of QSPRs which are often used in CAMD: (1) group contribution methods; (2) topological indices; and (3) signature descriptors. Next, in Section 3, we present the CAMD problem from a mathematical programming perspective, discussing various classes of the single-molecule design problem as well as CAMD problems considering mixtures of molecules and those involving the simultaneous design of a chemical product and the process it is a part of. In this section, several other important design considerations are presented, including a few important constraints to ensure practical solutions as well as the chemical feasibility of the designed structures. In Section 4, various solution techniques for the CAMD problem are discussed, including mathematical optimization strategies, decomposition methods, and heuristic approaches. Finally, in Section 5, we provide a diverse though non-exhaustive review of applications of CAMD problems.

Popular types of QSPRs in CAMD

The CAMD problem attempts to choose optimal (or simply good) molecules for some purpose from the space of theoretically possible chemical structures. At first glance, the CAMD problem must consider a very abstract chemical design space of atoms, bonds, aromaticity, structural isomers, electronic effects, etc. Though many of these features are certainly what gives molecules their specific properties and chemical functionality, they are difficult to build into any type of optimization scheme. This is primarily because there is no available immediate relationship between an arbitrary chemical structure and its performance or suitability regarding a specific application. In order to “rank” different structures and choose an optimal one, we must have some efficient way to quantify the properties and performance of each structure.

A second issue is the sheer size of the chemical search space. At the time of writing this article, the CAS registry [27] reports over 115 million unique organic and inorganic structures. This number only represents compounds which have been synthesized and cataloged, and it is already far too large for every structure to be considered in any type of trial-and-error design scheme. This number is also only a fraction of the theoretically possible chemical space, which some estimates indicate may contain more than $10^{60}$ unique molecules for small, drug-like structures [13]. Even with very efficient ways to estimate the performance of a certain structure, screening these structures using an enumeration strategy is far beyond current computational capacity. For this reason, we also need to relate the chemical space to a space that can be utilized for combinatorial optimization, allowing us to design over the massive search space far more efficiently.

CAMD practitioners have relied on semi-empirical quantitative structure property relationships (QSPRs) to address both of these issues. First, many semi-empirical methods delineate a clear connection between the abstract chemical space and a more practical space of quantitative properties. These methods are also often simple and can be applied to estimate properties very efficiently. Second, many of these methods break molecular structure into sub-molecular collections of atoms and bonds. These molecular sub-structures are assumed to dictate a molecule’s properties. Using these types of representations of the molecular space, combinatorial optimization can be directly applied to molecular design problems.

Group-contribution methods

The most commonly used QSPRs in CAMD are group contribution (GC) methods. These work under the assumption that a molecule’s properties can be predicted by the number of occurrences of various molecular sub-structures called “groups.” For example, we may think to represent the simple molecule propanol as a combination of the groups $-\text{CH}_{3}$, $-\text{CH}_{2}-$, and $-\text{OH}$. In this case, the dashes $(-)$ represent bonds to other groups. In its group representation, propanol would no longer be thought of as the connected alcohol molecule, but rather as some collection of its constituent groups. The group representation of propanol is shown in Fig. 1.

Being QSPRs, group contribution methods translate the group representation of a molecular structure into an estimate for some property $P$. To do this, group contribution methods define a vector $n$ that represents the number of occurrences of each of the groups. Assuming we only have the three groups shown above, propanol’s $n$ vector would look like $n=[1,2,1]$, where the entries in this vector represent the number of occurrences of the groups $-\text{CH}_{3}$, $-\text{CH}_{2}-$, and $-\text{OH}$, respectively. Each of these groups $g$ would also be associated with a coefficient $c_{g}$ which quantifies its affect or “contribution” to a particular property $P$. Properties are calculated as follows:

[TABLE]

The vector of coefficients $c$ comes from regression of the property $P$ over a large dataset of different molecules. To regress these parameters, the identity of all of the groups must be specified a priori. Returning to our example, one can easily imagine different sets of groups being used to describe propanol. For example, the groups $-\text{CH}_{3}$, $-\text{CH}_{2}-$, and $-\text{CH}_{2}\text{OH}$ also completely account for the atoms in propanol, and these may provide a better fit for the regression problem. For this reason, different group contribution methods to estimate different properties usually do not have consistent sets of groups. Finally, we note that the vector $n$ is generally much bigger as many group contribution methods contain 50-100 groups. Many group contribution methods make the additional assumption that groups cannot overlap, which means that $n$ is typically a sparse vector. A pictorial example of the usage of group contribution methods is given in Fig. 2. In this example, we apply a hypothetical GC method to a hypothetical molecule. We show how a molecular structure is decomposed into its constituent groups and provide a count of each of these groups. These counts constitute the elements of the vector $n$. The $n$ vector is paired with a hypothetical $c$ vector, and an example property is calculated.

One of the earliest examples of GC methods is from Benson and Buss [10], who are considered to be the originators of so-called “group increment theory.” Group increment theory, or Benson group increment theory (BGIT), is analogous to GC methods, but these terms may be more common in the physical chemistry literature. In the original 1958 paper [10], Benson and Buss proposed a simple group additivity scheme for the prediction of bond dissociation energies. Benson et al. [11] extended this work to account for a greater diversity of groups and to estimate heat capacities. Additional work from Cohen and Benson includes estimating heats of formation with group increment theory [31]. A large number of additional efforts have used Benson-like increments to estimate the same thermophysical properties, a very small sample of which are [36, 76, 144].

Another very popular GC method was devised by Joback and Reid [78]. This method extended the group increment idea to model many different properties with the same set of groups. The Joback and Reid model also included functional transformations for the original group increment summations. These altered the group contribution definition to the following:

[TABLE]

where $f$ represents some function of the inner product of the vectors $c$ and $n$. These functions $f$ appear in many group contribution methods and are important when predicted properties are not simple linear functions of the number of groups in a structure.

Perhaps the most widely-used GC method in CAMD is that of Marrero and Gani [113, 114]. This method, like the method of Constantinou and Gani [32] that predated it, provides another extension to the general form of GC methods in that it introduces multiple levels of groups to better capture proximity effects, meaning the effect of two or more groups which are close to one another in a molecular structure. As many GC methods are limited to groups of just a few atoms and bonds, many cannot differentiate between structures with different connectivity. The Marrero-Gani method, called the GC+ method, uses as a first-order approximation a normal group contribution method, where the groups belong to a set of primary groups $F$. An additional set of groups $S$ contains slightly larger sub-structures. Finally, a set of groups $T$ accounts for large groups and overarching molecular structural features. Unlike groups in the set $F$, groups in the sets $S$ and $T$ are allowed to overlap with each other. Using the hierarchical depiction of molecules with the GC+ method, a much clearer picture of a molecule is provided. The general form of the GC+ estimates is shown below:

[TABLE]

Group contribution methods have also often incorporated interaction terms [122, 123, 124, 125, 99, 140]. These are ways to include additional terms to account for the simultaneous presence of two (same or different) groups in a particular structure. For example, in predicting toxicity, one group in a molecule may lead to a simple metabolic pathway and therefore make many structures containing that group non-toxic. If that structure were also to have another group which is normally quite toxic, a GC method without interaction terms may not predict toxicity well. This would be because both groups would have an additive effect, meaning that there would be one highly non-toxic contribution and one highly toxic contribution. As a result, the molecule may be predicted to have an average toxicity when, in reality, the presence of the non-toxic group should outweigh the toxic group. Introducing an interaction term for these two groups accounts for the situation of their co-occurrence. In this example, this interaction term would likely remove whatever toxicity value was predicted by the toxic group. For more discussion of interaction terms in predicting toxicity, readers are referred to  [115]. Interaction terms generally take the form

[TABLE]

where $f_{I}$ usually represents multiplication, but can also represent other functions. Group contribution methods can also include some idea of structural features [122, 123, 124, 125, 113, 114]. These account for larger effects, typically at the molecule scale, such as aromatic ring substitution, cis/trans isomerism, aliphatic chain lengths, etc. These can typically be implemented as large groups but sometimes require special considerations. Table 1 provides several properties typically used in CAMD problems and a few GC methods to estimate them.

Strengths of GC methods. GC methods are useful in that they are very intuitive to use. They represent a chemical structure in terms of its functional components, very analogously to how chemists compare and analyze structures. GC methods are also able to easily represent a large and diverse chemical space as the groups can be combined in many different ways to produce a large variety of different structures. This is especially useful from a CAMD perspective. Finally, GC methods are easily translated into the mathematical formulations of CAMD problems as the inclusion and count of the groups (the vector $n$) are easily represented in the context of mathematical optimization.

Weaknesses of GC methods. There are a few shortcomings of modern GC methods. One is that many GC methods are unable to distinguish isomers from one another. As isomers can have very different properties, this represents a gap in the predictive power of GC methods. We note that some GC methods such as the GC+ methods are able to distinguish many isomers due to the inclusion of large groups. A second issue with GC methods is the lack of consistency in groups used to predict various properties. Though this has no major effect estimating these properties for a given structure, it becomes problematic for mathematical formulations of the CAMD problem. Finally, GC methods require specifying the set of groups prior to regressing the GC coefficients. Though many GC methods are quite accurate, there is no guarantee that the set of groups used best captures the property they model. Using different groups can sometimes drastically alter the predictive power of a GC model.

Topological indices

Chemical graph theory [15] is a field which became very influential in the 1970’s and has since been used to produce a large number of QSPRs. The basic idea of chemical graph theory is that the atoms and bonds which constitute a molecule can be thought of as nodes and edges in a graph. In general, we use $G=(V,E)$ to define a graph $G$, with vertex set $V$ and edge set $E$. Using this depiction of molecular structures, various properties of that graph, referred to here as topological indices, can be used as descriptors in QSPR models. More specifically, various topological indices are paired with regression coefficients and used to estimate properties in a similar way to GC methods.

Topological indices (TIs) can take many forms. They are defined as some function of the nodes and edges in a chemical graph, and one can easily see that there are a large number of possible functions even just considering standard molecular graph properties like degree counts for nodes, connectivity, atomic types, etc. One of the first topological indices used in chemical graph theory is the Wiener index [166]. The Wiener index attempts to describe the total distance between all atoms in the graph, as given by $d(v,v^{\prime})$, the graph theoretic distance between all pairs of vertices $v$ and $v^{\prime}$. The Wiener index $W(G)$ is defined as:

[TABLE]

While the Wiener index describes a graph in terms of its distances, another important consideration is how a graph is connected. To address this, an important class of topological indices called connectivity indices (CIs) was developed. Connectivity indices are widely used in CAMD and have been shown to be useful in QSPR applications [44]. The first connectivity indices were developed by Randić [142], who used these indices to account for the degree of branching in alkanes and to model enthalpy of fusion and vapor pressure. Randić defined an edge index to be:

[TABLE]

where $v$ and $v^{\prime}$ are two connected vertices in the chemical graph. This means that the atoms which correspond to $v$ and $v^{\prime}$ are connected by a chemical bond. Furthermore, $\delta_{v}$ and $\delta_{v^{\prime}}$ are the degrees of nodes $v$ and $v^{\prime}$. In the study of Randić, these degrees signified the number of $\sigma$ bonds a particular atom had to non-Hydrogen atoms (i.e., the number atomic neighbors in the hydrogen-suppressed graph), but they are sometimes defined differently for other connectivity indices. The connectivity index of the entire molecule (graph) was then given by

[TABLE]

where $E$ is the edge set of the graph. A calculation for the edge connectivity indices for a simple alkane is shown in Fig. 3.

${}^{1}\chi$ is the so-called first-order Randić connectivity index. An even simpler “connectivity index” exists which does not account for bonding at all. This is called the zeroth-order connectivity index, and is given below:

[TABLE]

Kier et al. [89], Hall et al. [65], and Murray et al. [121] were the first to apply these connectivity indices as descriptors in QSPR models. These models are typically linear in the descriptors, but many variations exist. An important step forward came from Kier et al. [90] who introduced higher-order connectivity indices. In general, these are defined for an index of order $i$ as:

[TABLE]

where $V_{C}^{i}$ is the set of all sets of $i$ connected vertices. Kier and Hall [87] developed additional connectivity indices to account for heteroatoms. These modified connectivity indices distinguished atoms by their valence electrons, leading to a new value for vertex degree, $\delta^{\text{V}}$, where the superscript indicates that valence is considered. For second period elements, Kier and Hall [87] define $\delta_{v}^{\text{V}}$ for a vertex $v$ as

[TABLE]

where $Z_{v}^{\text{V}}$ indicates the number of valence electrons for atom/vertex $v$ and $h_{v}$ is the number of hydrogens attached to $v$. For atoms in the third period and beyond, the main difference from this perspective is the number of core electrons. These are accounted for in the following:

[TABLE]

where $Z_{v}$ is the atomic number of atom $v$. These modified $\delta^{\text{V}}$ values define analogous connectivity indices ${}^{i}\chi^{\text{V}}$.

The connectivity index $\chi$ can be thought to primarily capture vertex adjacency and the local neighborhoods of every atom in a molecule. As such, it represents how atoms are connected and perhaps not the ensemble molecular structure and shape. Another topological index often used in CAMD that aims to address this is the so-called shape index, $\kappa$. The shape index accounts for features of the entire molecular structure as functions of the underlying graph architecture and counts of graph substructures. The main parameter in the calculation of $\kappa$ is the number of paths in the chemical graph of a certain length. We define ${}^{i}P$ to be the number of paths of length $i$ in a particular chemical graph.

Like other topological indices, the shape index maps this information onto a single value. There are three often used ${}^{i}\kappa$ values which are defined as:

[TABLE]

with ${}^{i}P_{\text{max}}$ and ${}^{i}P_{\text{min}}$ representing the maximum and minimum possible number of paths of length $i$ for a hypothetical molecule with an equivalent number of atoms. ${}^{i}P_{\text{max}}$ and ${}^{i}P_{\text{min}}$ can be easily derived from graph theoretic arguments. The formal expressions for these can be found in [63]. Several ${}^{i}\kappa$ values for an example molecule are calculated in Fig. 4.

In Table 2, we list a few references for QSPR models using topological indices. There are many, many examples of such models, and we note that this table is not meant to be exhaustive. For a more complete list of available topological-indices-based QSPRs, the reader is directed to the book of Devillers and Balaban [35].

Table 2 provides some idea of the diversity of applications of QSPR development with topological indices. We note that this table is inclined towards properties common in CAMD problems, so the properties listed may be more relevant to chemical engineering applications. However, we provide a few examples of the many applications of topological indices to modeling biological, environmental, and pharamacological properties. There have been many efforts to build models related to pharmaceutical properties [7, 45], so only a small subset is listed here. One additional interesting application of connectivity indices is in the prediction of coefficients for group contribution methods where the group is missing (the group is not in the descriptor space) [56].

Strengths of TIs. One of the main advantages of TIs is that they can discriminate between very similar structures, often in cases where GC methods cannot (e.g., isomers). This provides a more holistic picture of the molecule and can be very useful for certain design problems. For example, this may have potential applications for CAMD in areas where a structural feature of the to-be-designed compound is fixed a priori. Furthermore, since many TIs are a function of the entire graph, TIs reflect the entire nature of the molecular structure. This can have advantages over GC methods, which assume that each group provides a contribution independently of other groups in the structure (this is offset somewhat by GC interaction functions). Finally, TIs have been extensively applied to modeling pharmacological properties. The quality and volume of this literature means that TIs are very suitable to many pharmaceutically relevant CAMD problems.

Weaknesses of TIs. Though topological indices have been widely applied to QSPR, there are only limited examples in CAMD. Topological indices are usually not as generally applicable as GC methods, meaning that TI-based QSPRs are often restricted to a certain class of chemicals. For design purposes, this is problematic as it means that TI-based design problems can only consider that particular subset of the chemical search space. Furthermore, TIs represent graph-theoretic properties of the chemical graph, and many of these properties are not always readily understandable from a chemical perspective (although Randić and Zupan [143] have offered some interpretations of several TIs). Finally, TIs are more difficult to incorporate into CAMD than GC methods. We will discuss TI-based CAMD in a subsequent section, but such CAMD problems can sometimes face combinatorial difficulties and have only been demonstrated thus far on small design problems.

Signature descriptors

In a broad sense, GC methods capture the important subsets of atoms in a molecule while TIs rely on some function of the chemical graph. One QSPR method which has been shown to capture aspects of both GC- and TI-based methods is signature descriptors (SD). Signature descriptors are far younger than the other methods discussed, originating in 2002 from Visco et al. [163] and 2003 from Faulon et al. [47]. Like TIs, they conceive of chemical structures as the chemical graph. Rather than ascribe various values to a complete molecular graph, SDs retain all of the structural and connectivity information for every atom in a molecule.

Analogously to TIs, SDs define the chemical graph to be $G=(V,E)$. Standard SD methods also incorporate node coloring for every node $v$ via a coloring function $c_{G}(v)$ and the colors of each node $C_{v}$. This change is reflected in the slightly altered definition of the chemical graph as $G=(V,E,C,c_{G})$. These node colorings are intended to distinguish between different atoms as well as different types of the same atom. For example, it may be beneficial for the model to differentiate between an oxygen with two single bonds and an oxygen with a double bond. Additionally, one may want to distinguish hydrogens by what atom they are attached to, aromatic carbons from non-aromatic ones, atoms attached to aromatic rings, and many other chemical features. We note that the colorings of the nodes are one of the only subjective parts of SD models, and different coloring schemes can have significant effects on the performance of the models.

One important class of signature descriptors is known as atomic signatures. Given a certain atom in a chemical graph, its atomic signature represents all of the atoms that are within a certain distance, or height, from it. Varying the value of this distance gives rise to different atomic signatures. In the simplest case, with this distance set to 0, the atomic signature of an atom is simply that atom, colored in keeping with the coloring definition. More formally, for an atom (vertex) $v$, its atomic signature of height 0, ${}^{0}\alpha$ is given by:

[TABLE]

Of course, these atomic signatures of height 0 are not a robust set of descriptors. Considering higher values for height provides a more detailed picture of an atom’s environment. In general, for a height of $i$, the atomic signature of height $i\geq 1$ for a vertex $v$ is defined as:

[TABLE]

where $G_{S}$ defines a subgraph of $G$ that contains all vertices and bonds such that the distance between $v$ and any vertex in $G_{S}$ is at most $i$. Thus, an atomic signature of height one for a hypothetical atom $v$ defines $v$ and every atom bonded to $v$ as well as all connecting edges. A height two atomic signature defines $v$, every atom bonded to $v$, and every atom bonded to those atoms bonded to $v$ along with all necessary edges. Increasing the height of an atomic signature can thus be thought to add a layer of connected atoms. In graph-theoretic terms, increasing the $i$ value is equivalent to adding a layer to a breadth-first-search. A pictorial explanation is provided in Fig. 5. In this example, we assume that all atoms are colored by aromaticity and $sp$ hybridization and that hydrogens are colored by the atom they are attached to.

The atomic signatures can be thought of as a descriptor space for a molecular structure. A property, $P$, of a molecule can be estimated by all of its atomic signatures of up to a particular height. This QSPR has a familiar form:

[TABLE]

where $d$ is the index of the set of all atomic signatures and $D^{i}$ is the set of atomic SDs of height $i$. $c_{d}$ is a regression coefficient accounting for the “contribution” of each atomic signature to a certain property. ${}^{i}\alpha_{G}(d)$ represents the number of occurrences of atomic signature $d$. Using the atomic signatures, signatures of the entire molecule can also be generated.

A major advantage of signature descriptors is that they can be manipulated via simple functions to represent groups from GC methods as well as various TIs. This means that the large amount of QSPRs derived from GC methods and TI-based methods are accessible using SDs. For this reason, we do not provide a table of QSPRs using SDs because they can be—and often are—used to calculate properties via GC- and TI-based methods. A few examples of converting SDs to groups and TIs are given in [47].

Strengths of SDs. One of the main advantages of signature descriptors is that they have a small inherent bias as compared to GC methods and TIs. The only bias introduced in these descriptors is the choice of an atomic coloring scheme, or the choice of what defines a different type of the same atom. Furthermore, the inclusion of every atomic signature in SD-based QSPRs means there are no theoretical restrictions on the descriptor space defined for a molecular structure. There are also a variety of more modern signature descriptors capable of distinguishing stereoisomers. In the case of stereoisomers, certain SDs have far greater discriminative power than GC methods or TIs, which typically cannot differentiate stereoisomers. Furthermore, the equivalence of SDs to many TIs and groups for GC methods means that SDs can be used directly with these TI and GC QSPR models. This makes a large library of QSPR models accessible to SDs.

Weaknesses of SDs. Many QSPRs using SDs use all of the available atomic signatures of up to a certain height. This can quickly become an issue with the predictive power of SDs as models without sufficient training data may be overfit. A second concern with SDs is that the coloring scheme of the descriptors must be specified before the descriptors are used for QSPRs and CAMD problems. It is likely that some coloring schemes provide better results than others, and the best coloring scheme may not always be easy to determine. Finally, atomic signature descriptors always discriminate between identical atoms in different structural environments. This may again lead to issues of overfitting or not capturing the “true” chemical behavior of the system as sometimes it is better to model identical atoms with a general descriptor which is independent of the atom’s environment (although this can be captured to some degree with SDs of lower heights).

CAMD as an optimization problem

The various types of structure-property relationships discussed above can quickly and often accurately estimate properties from a structure. The application of QSPR techniques in this direction—predicting properties from structures—defines what is known as the “forward problem,” and is what QSPR techniques are generally intended for. CAMD can broadly be thought to consider the “reverse problem,” or the problem of predicting structures from properties. At first glance, there is no immediately obvious way of relating properties to a specific molecular structure. One issue is that there are so many structures to consider. A reasonable approach to the CAMD problem should be able to consider a large diversity of structures without running into significant computational difficulties. Another issue involves structural feasibility. Solutions to the CAMD problem must also be sensible molecular structures, meaning CAMD should produce structures that do not violate any inherent laws of chemical bonding. A final issue involves consistency. All of the QSPRs discussed (though non-overlapping GC methods are an exception) require that if certain features are present, so must be other features. For example, it would be unreasonable to design a structure that has four paths of length three but no paths of length two. In the case of GC methods or SDs, an example erroneous solution would contain one occurrence of the carboxylic acid (-C(=O)OH) group/signature and no occurrences of the carbonyl (-C(=O)-) group/signature. Note that this example assumes that overlap is allowed between these groups/signatures (although this is generally the case with signatures).

Mathematical optimization is the key to addressing all of these issues. Before introducing the problem in a general optimization formulation, we define a few important sets and variables. First, we assume that we have a vector of properties $p$ and a property value $p_{k}$ for each property $k$. The vector $n$ encapsulates relevant structural information of the designed molecules and is dependent on the type of QSPR chosen. In the case of GC methods, the value $n_{d}$ would represent the number of occurrences of each group $d$. For TIs, the $n$ vector may represent the number of topological features $d$ (edges, paths of certain lengths, etc.) from which TIs would be calculated. For SDs, this $n$ vector usually represents counts of various atomic signature descriptors $d$. The function $f$ then transforms this structural information into a property estimate using the appropriate QSPR relationship. In general, the CAMD problem can then be expressed as:

[TABLE]

In the above, (19) involves QSPR functions $f$ which estimate a vector of properties $p$ from attributes such as group counts, graph topological features, or atomic signatures. (20) (21) define general functions, $h_{1}$ and $h_{2}$, representing inequality and equality constraints on property values, desired structural features, process conditions, and a variety of other possibilities. For example, the presence of certain groups or structural features may necessitate changing a chemical process to accommodate these structures. These constraints can also eliminate certain groups/topological features/signatures from the solution space or require that they appear a certain number of times. The functions $h_{1}$ and $h_{2}$ can also account for design considerations such as system thermodynamics, cost, and a variety of system-specific interactions between $n$ and $p$. (22) and (23) define inequality and equality constraints which ensure structural feasibility. More specifically, the functions $s_{1}$ and $s_{2}$ determine if the vector $n$ is consistent with a molecular structure which can actually exist. These constraints prevent erroneous structures from being considered, eliminating compounds that violate atomic valences, are disjoint, have an incorrect number of aromatic atoms, etc. 24 and (25) set bounds on property values and $n$, respectively. Each property $k$ is bounded below by $p_{k}^{L}$ and above by $p_{k}^{U}$. Similarly, $n_{k}^{L}$ and $n_{k}^{U}$ define lower and upper bounds for $n_{k}$. Finally, (18) is a general objective function for the CAMD problem. $C(n,p)$ can define a number of possible functions. These functions somehow quantify the performance of a specific molecule based on its properties $p$ and perhaps its descriptors $n$. Of course, this $C$ function can either be minimized or maximized depending on the design problem.

Classes of the CAMD problem

Single molecule design

The single molecule design problem is the problem of determining a single, optimal structure for a particular purpose. For these problems, the structure of the compound is the only design consideration, meaning that the variables represented by $n$ above are the only degrees of freedom in the problem. Furthermore, it is assumed that there exists some ranking criterion with which to determine which structures are better than others. To make this section as general as possible, we assume that the ranking criterion can either be applied during the optimization procedure as the objective function or afterwards, evaluating the performance of each of a pool of candidate molecules. Though many classifications of these single molecule design problems are possible, we suggest three basic forms: (1) determining all feasible structures; (2) using an objective function which directly quantifies a molecule’s performance; and (3) designing a structure with properties as close as possible to certain property targets.

The feasibility problem. We begin by describing the case where there is no objective function used in the optimization problem. More formally, this is equivalent to solving the above CAMD formulation with the objective function equal to a constant. The solutions to this problem represent all molecules which satisfy the functions $h_{1}$ and $h_{2}$, are chemically sound structures, and do not violate the property or descriptor constraints. This is a useful type of problem to solve in CAMD when the performance function $C$ is not accessible or reliable. It may be the case that $C$ requires a complex simulation or experimental work. The performance function may also be inaccurate or the designer may not know exactly which properties to optimize. In these cases, these CAMD problems leverage the power of optimization to reduce the large number of possible compounds to a more manageable number. This smaller pool of compounds can then be investigated using high-order models or experiments.

Feasibility problems are sometimes difficult to distinguish from other types of problems. In many cases, a feasibility problem is solved at one of the beginning stages of the problem, and all of the feasible structures are then evaluated based on some ranking criteria. For these types of problems, there is no clear line between a feasibility problem and a problem with an objective function. Usually, the molecules designed by these problems are ultimately ranked, so many “feasibility problems” actually have an implicit objective function. A second complication is that many CAMD methodologies are designed to solve both feasibility problems and problems with objectives. We provide a few examples of feasibility problems here for GC methods [77], TIs [88, 91], and SDs [29]. We leave a more detailed discussion for the applications section of this article. A graphical example of the feasibility problem is given in Fig. 6.

Exact relationship between structures and performance. In many cases, it is possible to model a molecule’s performance directly as a function of its structure and properties. This is the most natural and common form of the CAMD objective because the performance function used to rank molecular structures is what is minimized or maximized. Using the exact performance function as the objective $C(n,p)$ guarantees that the solution to the CAMD problem is the optimal molecule for the application, at least as judged by the provided performance function. We also note that if the function $C(n,p)$ or any of the constraints is non-convex, the optimal molecule may only represent a local optimum. Using a global optimization algorithm will provide the best molecule for these non-convex problems.

An exact representation of the objective has inspired many numerical optimization strategies to solving the CAMD problem. Numerical optimization strategies are particularly advantageous in these cases because there is an exact algebraic relationship between the descriptor space $n$ and the performance function $C$. This enables the use of modern, state-of-the-art optimization techniques. Several of these will be discussed in an upcoming section.

There are many examples in the CAMD literature which fall into this category. Again, we will postpone the discussion of these topics until later sections of this document. We provide a small (and non-exhaustive) selection of relevant references in the meanwhile to offer some insight to the interested reader [127, 154, 146](GC), [20, 152] (TIs), and [28] (SDs). A graphical example of this type of CAMD problem is provided in Fig. 7.

Minimizing distance to property targets. Another class of single molecule design problems concerns finding a molecule with properties as close as possible to target values. These types of formulations have been used extensively in CAMD, often to design alternatives to molecules being used in practice. These molecules may need to be substituted for reasons of environmental friendliness, cost, or availability. Many other types of problems can be addressed when the ideal properties of a molecule are known. However, if these ideal properties do not represent a structure that is physically realistic, these approaches face some difficulty.

In general, these types of problems define the objective to be:

[TABLE]

where $p_{k}^{T}$ represents a targets for property $k$ and $w_{k}$ is a weight for the distance from from the estimated property value to its corresponding property target. The second norm shown above is the most common distance function in CAMD, but many other distance functions are possible.

We present a few demonstrative examples of CAMD with property target objective functions. Matsuda et al. [116] used groups to design ionic liquids based on conductivity and viscosity targets. Siddhaye et al. [153] addressed this problem with topological indices to design pharmaceutical products. Brown et al. [17] used signature descriptors to design polymers with specific properties. A graphical example of this type of CAMD problem is given in Fig. 8.

Mixture design

Real-world applications often demand a product with specifically tailored properties. This sometimes necessitates utilizing a mixture of compounds as no single compound possesses all of the necessary properties. Using CAMD techniques to simultaneously design two or more compounds for use in a blend/mixture is referred to as the mixture design problem. We note here that we alter a definition given in [4] for consistency with the prevailing definition of mixture design in the literature. All applications classified under mixture design can be assumed to design two or more structures simultaneously.

While the single-molecule design is difficult in many cases, the mixture design problem is much harder. The difficulties come from several sources: (1) the descriptor variables must now represent the descriptors of every unknown component in the mixture; (2) mixture properties must be calculated and included in the problem; (3) non-ideal mixture behavior must be considered in the form of thermodynamic relationships, activity coefficient models, or equations of state; (4) the design of mixtures also requires a determination of the amount of each component, so mole fractions must be considered. We provide a pictorial representation of the mixture design problem in Fig. 9. In this problem, we define the variable $x_{i}$ to represent the mole fraction of unknown component $i$. Furthermore, $q_{j}$ will represent the mixture property $j$. The formulation of the mixture design problem is similar to the single-molecule design problem. The objective function is altered to now include mixture properties:

[TABLE]

A few additional constraints are also necessary:

[TABLE]

In the above, (30) and (31) represent modified constraints from the previous formulation to include the mixture property variables $q$. Eq.(28) represents simple mixing functions or complex thermodynamic relationships which relate individual component property variables $p$, descriptor variables $n$, and mole fractions $x$ to the mixture property variables $q$. In many cases, some of these $q$ values represent activity coefficients. Eq. (29) ensures that mole fractions sum to 1. Finally, (32) places upper ($q_{j}^{U}$) and lower ($q_{j}^{L}$) bounds on mixture properties. We provide a few mixture design citations [54, 33, 19] and leave the rest for the applications section.

Integrated process and product design

Though many CAMD endeavors design products with the ultimate goal of being incorporated into an industrial process, few have explicitly considered the relationship between a particular structure and a process. This is especially important as process performance is typically very sensitive to the molecule(s) chosen. These problems are especially challenging from an optimization point of view because there is no easily discernible algebraic relationship between the descriptor variables $n$ and process variables. This problem also requires the introduction of process variables $\mu_{w}$, with $w$ defining the index over the process variables. To present the formulation for integrated process and product design, we modify the single-molecule design problem for simplicity.

First, the objective function now reflects process variables:

[TABLE]

A few additional constraints then account for the inclusion of process variables:

[TABLE]

In the above, (34) and (35) represent modified constraints from the previous formulation to include the process variables $\mu$. (36) is introduced to place upper ($\mu_{w}^{U}$) and lower ($\mu_{w}^{L}$) bounds on process variables. There is a similar formulation defined for process/product design problems which also consider mixtures. Again, we provide a small selection of references [41, 129, 82] and leave the remaining discussion for the applications section. One approach to the process design problem is illustrated in Fig. 10.

Common design features in the form of constraints

These constraints define any process criteria, design necessities, thermodynamic conditions, and any other features which may be important for a particular design problem. As these constraints reflect the diversity of applications of CAMD, we cannot go into detail about all of them here. Rather, we present a summary of design features and conditions which occur most commonly in CAMD.

A very common design feature in CAMD limits the number of groups which can occur in the solution. For example, constraints of this variety can place a limit on the number of molecular descriptors $d$:

[TABLE]

or a limit on the total number of descriptors:

[TABLE]

where $N^{U}$ represents the maximum number of structural features (groups, atoms, descriptors, topological features, etc.) in the designed molecule. Many design problems also set a lower bound on this summation, $N^{L}$, which is typically equal to 2. This ensures that more than one descriptor appear in the solution. Finally, some design problems seek solutions which are analogues within a given family of chemical structures. In this case, it may be necessary to include the descriptor variables corresponding to this structural family:

[TABLE]

where $D^{\text{FIX}}$ defines the set of descriptors which must occur in the solution and $n_{d}^{F}$ represents the number of descriptors $d$ which must occur for the structural family to be produced. This constraint is typically part of (25), where other descriptors may have lower bounds for other reasons.

Properties are also bounded similarly in equation (24). These constraints ensure that no process conditions, environmental regulations, toxicity thresholds, etc. are violated with the designed structures.

Mixture design problems or single component design problems involving mixtures typically have a number of common constraints. Most notably, these problems often require some idea of the activities of the chemical species in the mixture and thus necessitate incorporating an activity coefficient model. There are three such models often used in CAMD:

UNIFAC. UNIFAC [50] is a group contribution variant of the UNIQUAC equation [2] and has been used extensively in CAMD. UNIFAC is simple, accurate, and easy to incorporate into CAMD problems due to its use of groups. The original UNIFAC may lack accuracy outside of standard temperature and pressure ranges, but many extensions and re-parameterizations exist to address these cases. 2. 2.

SAFT. SAFT [25] is an accurate equation of state that is applicable in many temperature and pressure domains. It is gaining popularity in CAMD and several group contribution methods [159, 108, 109, 137, 133] have already been developed to estimate the SAFT parameters necessary for its use. 3. 3.

COSMO-RS and -SAC COSMO-RS [95, 97] and COSMO-SAC [107] are two post-processing methods for the COSMO solvation model [96]. Unlike other approaches, these COSMO-based models use electronic surface charge distributions that are calculated at the quantum chemistry level. These methods are now more amenable to CAMD due to the development of various COSMO-based group contribution methods [119, 120, 3].

Forms of the structural feasibility constraints

The constraints defined above in (22) and (23) ensure that the structural descriptors chosen with the variables $n_{d}$ are consistent with all other structural features such that they can be assembled to create a chemically feasible molecule. The three different types of QSPRs given above require different sets of these constraints, so this discussion will be divided into three parts, beginning with group contribution methods.

Connectivity constraints for GC methods

If using GC methods to solve CAMD problems, the variable $n_{d}$ typically represents the number of occurrences of the group $d$ in the solution. Many structural constraints in CAMD also consider the valence of each group, $\Phi_{d}$, which is simply the number of bonds a group requires to satisfy its valence electron requirements. For example, the group $-\text{CH}_{\text{1}}\text{(Cl)}-$ requires 2 external bonds ($\Phi_{\text{CH}_{\text{1}}\text{(Cl)}}=2$) as two of carbon’s four required bonds are accounted for by hydrogen and chlorine atoms. One of the most widely-used rules for structural feasibility using valences comes from Odele and Macchietto [127]:

[TABLE]

where m is defined by:

[TABLE]

These constraints alone, while appropriate for most molecules, still allow for certain sets of groups which cannot be joined to form feasible molecules. For example, a solution of 1 $>\text{CH}_{\text{0}}<$ ($\Phi_{\text{CH}_{\text{0}}}=4$) group and 2 -Br ($\Phi_{\text{Br}}=1$) groups is feasible for the above constraints with $m=0$. To address these situations, Odele and Macchietto also define a constraint to ensure that there are enough groups to meet the valence requirements of every group. This means that to include a group in the solution with a valence requirement of 4, there must also be at least 4 additional groups. This is captured in

[TABLE]

where $d^{\prime}$ is also an index over descriptors. These constraints can be generalized with integer variables $N_{R}^{\text{arom}}$ and $N_{R}^{\text{ali}}$, which represent the number of aromatic rings and the number of aliphatic rings, respectively. Furthermore, we introduce $G_{ali}$ and $G_{arom}$ to be the sets of groups which have any open valences which are aliphatic and aromatic, respectively. We define these valence counts to be $\Phi^{ali}$ for aliphatics and $\Phi^{arom}$ for aromatics. Note that the same group can appear in both sets. We also define the parameter $A^{\text{arom}}_{d}$ to represent the number of aromatic atoms in a group $d$ and the parameter $\rho_{d}$ to represent the number of available aliphatic attachment points for every aromatic atom in group $d$. Below we provide a slightly more general formulation considering aromatic rings:

[TABLE]

We note that in this example, the aromatic rings are assumed to be benzylic and not attached to another aromatic ring (biphenyls and fused aromatics). An example of these basic structural constraints for GC methods are given in Fig. 11. A few simple extensions allow these constraints to account for all aromatic rings. These modified Odele-Macchietto constraints work in cases where all groups are bonded to each other with only single or aromatic bonds. More complex connectivity constraints, accounting for cases where groups may be double or triple bonded to one another are given in [146].

Connectivity constraints for TI methods