Analysis of the activities and properties of chemical compounds is important not only for chemical science but also for biological science because chemical compounds play important roles in metabolic and many other pathways. Computational prediction of chemical activities from their structural data has been studied for several decades under the name of quantitative structure activity relationship (QSAR) [1]. In addition to traditional regression-based methods, various machine learning methods have been applied to QSAR [2, 3]. Recently, neural networks and deep-learning technologies have extensively been applied to QSAR [4].

Inference of chemical structures with desired chemical activities under some constraints is also important because of its potential applications to drug design, and the problem has been studied under the name of inverse quantitative structure activity relationship (inverse QSAR). Chemical compounds are commonly represented by undirected graphs called chemical graphs in which vertices and edges correspond to atoms and chemical bonds, respectively. Due to the difficulty of directly handling chemical graphs in both QSAR and inverse QSAR, chemical compounds are usually represented as vectors of integer or real numbers, which are called descriptors in chemoinformatics and correspond to feature vectors in machine learning. In inverse QSAR, one major approach is to first infer feature vectors from given chemical activities and constraints, and then reconstruct chemical structures from these feature vectors [5, 6, 7]. However, the reconstruction itself is not an easy task because the number of possible chemical graphs is huge. For example, the number of chemical graphs with up to 30 atoms (vertices) C, N, O, and S may exceed ${10}^{60}$ [8]. Indeed, the problem to infer a chemical graph from a given feature vector is known as a computationally difficult problem (precisely, NP-hard) except for some simple cases [9]. Most existing methods for inverse QSAR do not guarantee exact or optimal solutions due to these inherent difficulties.

Recently, artificial neural networks (ANNs), in particular, graph convolutional networks [10] are extensively used for inverse QSAR. For example, recurrent neural networks [11, 12], variational autoencoders [13], grammar variational autoencoders [14], invertible flow models [15, 16], and generative adversarial networks [17] have been applied. However, these methods do not yet guarantee exact or optimal solutions.

Akutsu and Nagamochi [18] proved that the computation process of a given ANN can be simulated as a mixed integer linear programming (MILP). Based on this result, a novel framework for inferring a set of chemical graphs has been developed [19, 20], which is illustrated in Fig. 1. This framework consists of two phases. In the first phase, it constructs a prediction function and in the second phase, it infers a chemical graph. There are three stages in the first phase of the framework. In Stage 1, a chemical property $\pi$ and a class $𝒢$ of graphs are selected, and a property function $a$ is defined so that $a(ℂ)$ is the value of $\pi$ for a compound $ℂ\in 𝒢$. Then we collect a data set $D_{\pi }$ of chemical graphs in $𝒢$ such that $a(ℂ)$ is available for every $ℂ\in D_{\pi }$. In Stage 2, a feature function $f:𝒢\to {ℝ}^K$ for a positive integer $K$ is introduced. In Stage 3, a prediction function $\eta$ is constructed with an ANN $𝒩$ that, given a vector $x\in {ℝ}^K$, returns a value $y=\eta (x)\in ℝ$ so that $\eta (f(ℂ))$ serves as a predicted value to $a(ℂ)$ of $\pi$ for each $ℂ\in D_{\pi }$. Given a target chemical value $y^{*}$, the second phase consists the next two phases to infer chemical graphs ${ℂ}^{*}$ with $\eta \left(f\left({ℂ}^{*}\right)\right)=y^{*}$. A feature function $f$ and a prediction function $\eta$ are obtained in the first phase, and we call an additional constraint on the substructures of target chemical graphs a topological specification. In Stage 4, the following two MILP formulations are designed:

- MILP $ℳ(x,y;{𝒞}_1)$ with a set ${𝒞}_1$ of linear constraints on variables $x$ and $y$ (and some other auxiliary variables) simulates the process of computing $y:=\eta (x)$ from a vector $x$; and

- MILP $ℳ(g,x;{𝒞}_2)$ with a set ${𝒞}_2$ of linear constraints on variable $x$ and a variable vector $g$ that represents a chemical graph $ℂ$ (and some other auxiliary variables) simulates the process of computing $x:=f(ℂ)$ from a chemical graph $ℂ$and chooses a chemical graph $ℂ$ that satisfies the given topological specification $\sigma$.

Given a target value $y^{*}\in ℝ$, the combined MILP $ℳ(g,x,y;{𝒞}_1,{𝒞}_2)$ is solved to find a feature vector $x^{*}\in {ℝ}^K$ and a chemical graph ${ℂ}^{†}$ that satisfies the specification $\sigma$ such that $f({ℂ}^{†})=x^{*}$ and $\eta (x^{*})=y^{*}$ (where if the MILP is infeasible then this suggests that such a desired chemical graph does not exist). In Stage 5, by using the inferred chemical graph ${ℂ}^{†}$, we generate other chemical graphs ${ℂ}^{*}$ such that $\eta (f({ℂ}^{*}))=y^{*}$.

Stage 4 MILP formulations to infer chemical graphs with cycle index 0, 1 and 2 are proposed in [20, 21, 22, 23], respectively, but no sophisticated topological specification was available yet. A restricted class of acyclic graphs that is characterized by an integer $\rho$, called a “branch-parameter” is introduced by Azam et al. [21]. This restricted class still covers most of the acyclic chemical compounds in the database. Akutsu and Nagamochi [24] extended the idea to define a restricted class of cyclic graphs, called “$\rho$-lean cyclic graphs” and introduced a set of flexible rules for describing a topological specification. Tanaka et al. [25] used a decision tree instead of ANNs to construct a prediction function $\eta$ in Stage 3 in the framework and an MILP $ℳ(x,y;{𝒞}_1)$ that simulates the computation process of a decision tree.

Recently Shi et al. [26] proposed a new model to deal with an arbitrary graph in the framework called a two-layered model to represent the feature of a chemical graph. Also, the set of rules for describing a topological specification in [27] was refined so that a prescribed structure can be included in both of the acyclic and cyclic parts of a chemical graph $ℂ$. In this model, a chemical graph $ℂ$ with an integer $\rho \ge 1$, we consider two parts, namely, the exterior and the interior of the hydrogen-suppressed graph $⟨ℂ⟩$ that is obtained by removing hydrogen from $ℂ$. The exterior consists of maximal acyclic induced subgraphs with height at most $\rho$ in $⟨ℂ⟩$ and the interior is the connected subgraph of $⟨ℂ⟩$ obtained by excluding the exterior. Shi et al. [26] also defined a feature vector $f(ℂ)$ of a chemical graph $ℂ$ as a combination of the frequency of adjacent atom pairs in the interior and the frequency of chemical acyclic graphs among the set of chemical rooted trees $T_u$ rooted at interior-vertices $u$. Recently, Tanaka et al. [25] extended the model in order to directly treat a chemical graph with hydrogens so that the feature of the exterior can be represented with more variety of chemical rooted trees.

The contribution of this paper is as follows. Firstly, we make a slight modification to a model of chemical graphs proposed by Tanaka et al. [25] so that we can treat a chemical element with multi-valence such as sulfur S and a chemical graph with cations and anions. Then, we consider the prediction function. One of the most important factors in the framework is the quality of a prediction function $\eta$ constructed in Stage 3. Also, overfitting is pointed out as a major issue in ANN-based approaches for QSAR because many parameters need to be optimized for ANNs [4]. In this paper, to construct a prediction function in Stage 3, we use linear regression instead of ANNs or decision trees. A learning algorithm for an ANN may not find a set of weights and biases that minimizes the error function since the algorithm simply iterates modification of the current weights and biases until it terminates at a local minimum value, and linear regression is much simpler than ANNs and decision trees and thereby we regard the performance of a prediction function by linear regression as the basis for other more sophisticated machine learning methods. In this paper, we derive an MILP formulation $ℳ(x,y;{𝒞}_1)$ in Stage 4 to simulate the computation process of a prediction function by linear regression. For an MILP formulation $ℳ(g,x;{𝒞}_2)$ that represents a feature function $f$ and a specification $\sigma$ in Stage 4, we can use the same formulation proposed by Tanaka et al. [25] with a slight modification (the detail of the MILP $ℳ(g,x;{𝒞}_2)$ can be found in Supplementary Material). In Stage 5, we can also use the dynamic programming algorithm due to Tanaka et al. [25] with a slight modification to generate target chemical graphs ${ℂ}^{*}$ and the details are omitted in this paper.

We implemented the framework based on the refined two-layered model and a prediction function by linear regression. The results of our computational experiments reveal a set of chemical properties to which a prediction function constructed by linear regression on our feature function performs well, in comparison with ANNs in our previous work. We also observe that chemical graphs with up to 50 non-hydrogen atoms can be inferred by the proposed method.

The paper is organized as follows. Section 2 introduces some notions and terminologies on graphs, modeling of chemical compounds and our choice of descriptors. Section 3 describes our modification to the two-layered model. Section 4 reviews the idea of linear regression and formulates an MILP $ℳ(x,y;{𝒞}_1)$ that simulates the computing process of a prediction function constructed by linear regression. Section 5 reports the results of some computational experiments conducted for 18 chemical properties such as vapor density and optical rotation. Section 6 gives conclusions with future work. Some technical details are given in Supplementary Material.

In this section, we review some notions and terminologies related to graphs, modeling of chemical compounds introduced by Tanaka et al. [25] and our choice of descriptors.

Let $ℝ$, ${ℝ}_{+}$, $\mathbb{Z}$ and ${\mathbb{Z}}_{+}$ denote the sets of reals, non-negative reals, integers and non-negative integers, respectively. For two integers $a$ and $b$, let $[a,b]$ denote the set of integers $i$ with $a\le i\le b$.

Graph Given a graph $G$, let $V(G)$ and $E(G)$ denote the sets of vertices and edges, respectively. For a subset $V^{\prime }\subseteq V(G)$ (resp., $E^{\prime }\subseteq E(G))$ of a graph $G$, let $G-V^{\prime }$ (resp., $G-E^{\prime }$) denote the graph obtained from $G$ by removing the vertices in $V^{\prime }$ (resp., the edges in $E^{\prime }$), where we remove all edges incident to a vertex in $V^{\prime }$ in $G-V^{\prime }$. An edge subset $E^{\prime }\subseteq E(G)$ in a connected graph $G$ is called separating (resp., non-separating) if $G-E^{\prime }$ becomes disconnected (resp., $G-E^{\prime }$ remains connected). The rank$\mathrm{r}(G)$ of a connected graph $G$ is defined to be the minimum $|F|$ of an edge subset $F\subseteq E(G)$ such that $G-F$ contains no cycle, where $\mathrm{r}(G)=|E(G)|-|V(G)|+1$. Observe that $\mathrm{r}(G-E^{\prime })=\mathrm{r}(G)-|E^{\prime }|$ holds for any non-separating edge subset $E^{\prime }\subseteq E(G)$. An edge $e=u_1{u}_2\in E(G)$ in a connected graph $G$ is called a bridge if $\{e\}$ is separating, i.e., $G-e$ consists of two connected graphs $G_i$ containing vertex $u_i$, $i=1,2$. For a connected cyclic graph $G$, an edge $e$ is called a core-edge if it is in a cycle of $G$ or is a bridge $e=u_1{u}_2$ such that each of the connected graphs $G_i$, $i=1,2$, of $G-e$ contains a cycle. A vertex incident to a core-edge is called a core-vertex of $G$. A path with two end-vertices $u$ and $v$ is called a $u,v$-path.

A vertex designated in a graph $G$ is called a root. In this paper, we designate at most two vertices as roots, and denote by $\mathrm{Rt}(G)$ the set of roots of $G$. We call a graph $G$rooted (resp., bi-rooted) if $|\mathrm{Rt}(G)|=1$ (resp., $|\mathrm{Rt}(G)|=2$), where we call $G$unrooted if $\mathrm{Rt}(G)=\mathrm{\varnothing }$.

For a graph $G$ possibly with roots, a leaf-vertex is defined to be a non-root vertex $v\in V(G)\backslash \mathrm{Rt}(G)$ with degree 1. Call the edge $uv$ incident to a leaf vertex $v$ a leaf-edge, and denote by $V_{\text{leaf}}(G)$ and $E_{\text{leaf}}(G)$ the sets of leaf-vertices and leaf-edges in $G$, respectively. For a graph or a rooted graph $G$, we define graphs $G_i,i\in {\mathbb{Z}}_{+}$ obtained from $G$ by removing the set of leaf-vertices $i$ times so that

(1)$$G_0:=G;G_{i+1}:=G_i-V_{\text{leaf}}\left(G_i\right),$$

where we call a vertex $v\in V_{\text{leaf}}\left(G_k\right)$ a leaf $k$-branch and we say that a vertex $v\in V_{\text{leaf}}\left(G_k\right)$ has $\text{height ht}(v)=k$ in $G$. The $\text{height ht}(T)$ of a rooted tree $T$ is defined to be the maximum of $\mathrm{ht}(v)$ of a vertex $v\in V(T)$. For an integer $k\ge 0$, we call a rooted tree $T$ $k$-lean if $T$ has at most one leaf $k$-branch. For an unrooted cyclic graph $G$, we regard that the set of non-core-edges in $G$ induces a collection $𝒯$ of trees each of which is rooted at a core-vertex, where we call $G$ $k$-lean if each of the rooted trees in $𝒯$ is $k$-lean.

This section reviews the idea of the two-layered model introduced by Shi et al. [26], and describes our modifications to the model.

Let $ℂ=(H,\alpha ,\beta )$ be a chemical graph and $\rho \ge 1$ be an integer, which is called a branch-parameter.

A two-layered model of $ℂ$ introduced by Shi et al. [26] is a partition of the hydrogen-suppressed chemical graph $⟨ℂ⟩$ into an “interior” and an “exterior” in the following way. We call a vertex $v\in V(⟨ℂ⟩)$ (resp., an edge $e\in E(⟨ℂ⟩))$ of $G$ an exterior-vertex (resp., exterior-edge) if (resp., $e$ is incident to an exterior-vertex) and denote the sets of exterior-vertices and exterior-edges by $V^{\mathrm{ex}}(ℂ)$ and $E^{\mathrm{ex}}(ℂ)$, respectively and denote $V^{\mathrm{int}}(ℂ)=V(⟨ℂ⟩)\backslash V^{\mathrm{ex}}(ℂ)$ and $E^{\mathrm{int}}(ℂ)=E(⟨ℂ⟩)\backslash E^{\mathrm{ex}}(ℂ)$, respectively. We call a vertex in $V^{\mathrm{int}}(ℂ)$ (resp., an edge in $E^{\mathrm{int}}(ℂ)$) an interior-vertex (resp., interior-edge). The set $E^{\mathrm{ex}}(ℂ)$ of exterior-edges forms a collection of connected graphs each of which is regarded as a rooted tree $T$ rooted at the vertex $v\in V(T)$ with the maximum $\mathrm{ht}(v)$. Let ${𝒯}^{\mathrm{ex}}(⟨ℂ⟩)$ denote the set of these chemical rooted trees in $⟨ℂ⟩$. The interior${ℂ}^{\text{int}}$ of $ℂ$ is defined to be the subgraph $(V^{\mathrm{int}}(ℂ),E^{\mathrm{int}}(ℂ))$ of $⟨ℂ⟩$.

Fig. 2 illustrates an example of a hydrogen-suppressed chemical graph $⟨ℂ⟩$. For a branch-parameter $\rho =2$, the interior of the chemical graph $⟨ℂ⟩$ in Fig. 2 is obtained by removing the set of vertices with degree 1 $\rho =2$ times; i.e., first remove the set $V_1=\{w_1,w_2,\mathrm{\dots },w_{14}\}$ of vertices of degree 1 in $⟨ℂ⟩$ and then remove the set $V_2=\{w_{15},w_{16},\mathrm{\dots },w_{19}\}$ of vertices of degree 1 in $⟨ℂ⟩-V_1$, where the removed vertices become the exterior-vertices of $⟨ℂ⟩$.

For each interior-vertex $u\in V^{\mathrm{int}}(ℂ)$, let $T_u\in {𝒯}^{\mathrm{ex}}(⟨ℂ⟩)$ denote the chemical tree rooted at $u$ (where possibly $T_u$ consists of vertex $u$) and define the $\rho$-fringe-tree$ℂ[u]$ to be the chemical rooted tree obtained from $T_u$ by putting back the hydrogens originally attached to $T_u$ in $ℂ$. Let $𝒯(ℂ)$ denote the set of $\rho$-fringe-trees $ℂ[u],u\in V^{\mathrm{int}}(ℂ)$. Fig. 3 illustrates the set $𝒯(ℂ)=\{ℂ\left[u_i\right]\mid i\in [1,28]\}$ of the 2-fringe-trees of the example $ℂ$ in Fig. 2.

Feature Function We extend the feature function of a chemical graph $ℂ$ introduced by Tanaka et al. [25]. The feature of an interior-edge $e=uv\in E^{\mathrm{int}}(ℂ)$ such that $\alpha (u)=\mathrm{a}$, ${\mathrm{deg}}_{⟨ℂ⟩}(u)=d$, $\alpha (v)=\mathrm{b}$, ${\mathrm{deg}}_{⟨ℂ⟩}(v)=d^{\prime }$ and $\beta (e)=m$ is represented by a tuple $(\mathrm{a}d,\mathrm{b}{d}^{\prime },m)$, which is called the edge-configuration of the edge $e$, where we call the tuple $(\mathrm{a},\mathrm{b},m)$ the adjacency-configuration of the edge $e$.

For an integer $K$, a feature vector $f(ℂ)$ of a chemical graph $ℂ$ is defined by a feature function$f$ that consists of $K$ descriptors. We call ${ℝ}^K$ the feature space.

Tanaka et al. [25] defined a feature vector $f(ℂ)\in {ℝ}^K$ to be a combination of the frequency of edge-configurations of the interior-edges and the frequency of chemical rooted trees among the set of chemical rooted trees $ℂ[u]$ over all interior-vertices $u$. In this paper, we introduce the rank and the adjacency-configuration of leaf-edges as new descriptors in a feature vector of a chemical graph. See Supplementary Material for a full description of descriptors used in Stage 2 of the framework.

For an integer $p\ge 1$ and a vector $x\in {ℝ}^p$, the $j$-th entry of $x$ is denoted by $x(j),j\in [1,p]$.

Let $D$ be a data set of chemical graphs $ℂ$ with an observed value $a(ℂ)\in ℝ$, where we denote by $a_i=a({ℂ}_i)$ for an indexed graph ${ℂ}_i$.

Let $f$ be a feature function that maps a chemical graph $ℂ$ to a vector $f(ℂ)\in {ℝ}^K$ where we denote by $x_i=f({ℂ}_i)$ for an indexed graph ${ℂ}_i$. For a prediction function $\eta :{ℝ}^K\to ℝ$, define an error function

(4)$$\mathrm{Err}(\eta ;\mathrm{D})\triangleq \sum _{{ℂ}_{\mathrm{i}}\in \mathrm{D}}{\left({\mathrm{a}}_{\mathrm{i}}-\eta \left(\mathrm{f}\left({ℂ}_{\mathrm{i}}\right)\right)\right)}^2=\sum _{{ℂ}_{\mathrm{i}}\in \mathrm{D}}{\left({\mathrm{a}}_{\mathrm{i}}-\eta \left({\mathrm{x}}_{\mathrm{i}}\right)\right)}^2,$$

and define the coefficient of determination $R^2(\eta ,D)$ to be

(5)$${\mathrm{R}}^2(\eta ,\mathrm{D})\triangleq 1-\frac{\mathrm{Err}(\eta ;\mathrm{D})}{{\sum }_{{ℂ}_{\mathrm{i}}\in \mathrm{D}}{\left({\mathrm{a}}_{\mathrm{i}}-\tilde{\mathrm{a}}\right)}^2}\text{for}\tilde{\mathrm{a}}=\frac{1}{|\mathrm{D}|}\sum _{ℂ\in \mathrm{D}}\mathrm{a}(ℂ).$$

For a feature space ${ℝ}^K$, a hyperplane is denoted by a pair $(w,b)$ of a vector $w\in {ℝ}^K$ and a real $b\in ℝ$. Given a hyperplane $(w,b)\in {ℝ}^{K}w$, a prediction function ${\eta }_{w,b}:{ℝ}^K\to ℝ$ is defined by setting

(6)$${\eta }_{w,b}(x)\triangleq w\cdot x+b=\sum _{j\in [1,K]}w(j)x(j)+b.$$

We wish to find a hyperplane $(w,b)$ that minimizes the error function $\mathrm{Err}({\eta }_{w,b};D)$. In many cases, a feature vector $f$ contains descriptors that do not play an essential role in constructing a good prediction function. When we solve the minimization problem, the entries $w(j)$ for some descriptors $j\in [1,K]$ in the resulting hyperplane $(w,b)$ become zero, which means that these descriptors were not necessarily important for finding a prediction function ${\eta }_{w,b}$. It is proposed that solving the minimization with an additional penalty term $\tau$ to the error function often results in more number of entries $w(j)=0$, reducing a set of descriptors necessary for defining a prediction function ${\eta }_{w,b}$. For an error function with such a penalty term, a Ridge function $\frac{1}{2|D|}\mathrm{Err}({\eta }_{w,b};D)+\lambda \left[{\sum }_{j\in [1,K]}w{(j)}^2+b^2\right]$ [28] and a Lasso function $\frac{1}{2|D|}\mathrm{Err}({\eta }_{w,b};D)+\lambda \left[{\sum }_{j\in [1,K]}|w(j)|+|b|\right]$ [29] are known, where $\lambda \in ℝ$ is a given real number.

Given a prediction function ${\eta }_{w,b}$, we can simulate a process of computing the output ${\eta }_{w,b}(x)$ for an input $x\in {ℝ}^K$ as an MILP $ℳ(x,y;{𝒞}_1)$ in the framework. By solving such an MILP for a specified target value $y^{*}$, we can find a vector $x^{*}\in {ℝ}^K$ such that ${\eta }_{w,b}(x^{*})=y^{*}$. Instead of specifying a single target value $y^{*}$, we use lower and upper bounds ${\underset{¯}{y}}^{*},{\overline{y}}^{*}\in ℝ$ on the value $a(ℂ)$ of a chemical graph $ℂ$ to be inferred. We can control the range between ${\underset{¯}{y}}^{*}$ and ${\overline{y}}^{*}$ for searching a chemical graph $ℂ$ by setting ${\underset{¯}{y}}^{*}$ and ${\overline{y}}^{*}$ to be close or different values. A desired MILP is formulated as follows.

$ℳ(x,y;{𝒞}_1)$: An MILP formulation for the inverse problem to prediction function.

constants:

- A hyperplane $(w,b)$ with $w\in {ℝ}^K$ and $b\in ℝ$;

- Real values ${\underset{¯}{y}}^{*},{\overline{y}}^{*}\in ℝ$ such that ;

- A set $I_{\mathbb{Z}}$ of indices $j\in [1,K]$ such that the $j$-th descriptor ${\mathrm{dcp}}_j(ℂ)$ is always an integer;

- A set $I_{+}$ of indices $j\in [1,K]$ such that the $j$-th descriptor ${\mathrm{dcp}}_j(ℂ)$ is always non-negative;

- $\mathrm{\ell }(j),u(j)\in ℝ$, $j\in [1,K]$: lower and upper bounds on the $j$-th descriptor;

variables:

- Non-negative integer variable $x(j)\in {\mathbb{Z}}_{+}$, $j\in I_{\mathbb{Z}}\cap I_{+}$;

- Integer variable $x(j)\in \mathbb{Z},j\in I_{\mathbb{Z}}\backslash I_{+}$;

- Non-negative real variable $x(j)\in {ℝ}_{+},j\in I_{+}\backslash I_{\mathbb{Z}}$;

- Real variable $x(j)\in ℝ,j\in [1,K]\backslash (I_{\mathbb{Z}}\cup I_{+})$;

constraints:

(7)$\mathrm{\ell }(j)\le x(j)\le u(j),j\in [1,K];$
${\underset{¯}{y}}^{*}\le {\displaystyle \sum _{j\in [1,K]}}w(j)x(j)+b\le {\overline{y}}^{*}$

objective function: none.

The number of variables and constraints in the above MILP formulation is $O(K)$. It is not difficult to see that the above MILP is an NP-hard problem. The entire MILP for Stage 4 consists of the two MILPs $ℳ(x,y;{𝒞}_1)$ and $ℳ(g,x;{𝒞}_2)$ with no objective function. The latter represents the computation process of our feature function $f$ and a given topological specification. See Supplementary Material for the details of MILP $ℳ(g,x;{𝒞}_2)$.

We implemented our method of Stages 1 to 5 for inferring chemical graphs under a given topological specification and conducted experiments to evaluate the computational efficiency. We executed the experiments on a PC with Processor: Core i7-9700 (3.0 GHz; 4.7 GHz at the maximum) and Memory: 16 GB RAM DDR4.

Results on Phase 1. We have conducted experiments of linear regression for 37 chemical properties among which we report the following 18 properties to which the test coefficient of determination ${\mathrm{R}}^2$ attains at least 0.8: octanol/water partition coefficient (Kow), heat of combustion (Hc), vapor density (Vd), optical rotation (OptR), electron density on the most positive atom (EDPA), melting point (Mp), heat of atomization (Ha), heat of formation (Hf), internal energy at 0K (U0), energy of lowest unoccupied molecular orbital (Lumo), isotropic polarizability (Alpha), heat capacity at 298.15K (Cv), solubility (Sl), surface tension (SfT), viscosity (Vis), isobaric heat capacities in liquid phase (IhcLiq), isobaric heat capacities in solid phase (IhcSol) and lipophilicity (Lp).

We used data sets provided by HSDB from PubChem [30]for Kow, Hc, Vd and OptR, M. Jalali-Heravi and M. Fatemi [31] for EDPA, Roy and Saha [32] for Mp, Ha and Hf, MoleculeNet [33] for U0, Lumo, Alpha, Cv and Sl, Goussard et al. [34] for SfT and Vis, R. Naef [35] for IhcLiq and IhcSol, and Figshare [36] for Lp.

Properties U0, Lumo, Alpha and Cv share a common original data set $D^{*}$ with more than 130,000 compounds, and we used a set $D_{\pi }$ of 1,000 graphs randomly selected from $D^{*}$ as a common data set of these four properties $\pi$ in this experiment.

Stages 1, 2 and 3 in Phase 1 are implemented as follows.

Stage 1. We set a graph class $𝒢$ to be the set of all chemical graphs with any graph structure, and set a branch-parameter $\rho$ to be 2.

For each of the properties, we first select a set $\mathrm{\Lambda }$ of chemical elements and then collect a data set $D_{\pi }$ on chemical graphs over the set $\mathrm{\Lambda }$ of chemical elements. During construction of the data set $D_{\pi }$, chemical compounds that do not satisfy one of the following are eliminated: the graph is connected, the number of non-hydrogen neighbors of each atom is at most four, and the number of carbon atoms is at least four.

Table 1 shows the size and range of data sets that we prepared for each chemical property in Stage 1, where we denote the following:

- $\mathrm{\Lambda }$: the set of elements used in the data set $D_{\pi }$; $\mathrm{\Lambda }$ is one of the following 11 sets: ${\mathrm{\Lambda }}_1=\{\mathrm{H},\mathrm{C},\mathrm{O}\}$; ${\mathrm{\Lambda }}_2=\{\mathrm{H},\mathrm{C},\mathrm{O},\mathrm{N}\}$; ${\mathrm{\Lambda }}_3=\{\mathrm{H},\mathrm{C},\mathrm{O},{\mathrm{S}}_{(2)}\}$; ${\mathrm{\Lambda }}_4=\{\mathrm{H},\mathrm{C},\mathrm{O},\mathrm{S}\text{i}\}$; ${\mathrm{\Lambda }}_5=\{\mathrm{H},\mathrm{C},\mathrm{O},\mathrm{N},\mathrm{C}\text{l},{\mathrm{P}}_{(3)},{\mathrm{P}}_{(5)}\}$; ${\mathrm{\Lambda }}_6=\{\mathrm{H},\mathrm{C},\mathrm{O},\mathrm{N},{\mathrm{S}}_{(2)},\mathrm{F}\}$; ${\mathrm{\Lambda }}_7=\{\mathrm{H},\mathrm{C},\mathrm{O},\mathrm{N},{\mathrm{S}}_{(2)},{\mathrm{S}}_{(6)},\mathrm{C}\text{l}\}$; ${\mathrm{\Lambda }}_8=\{\mathrm{H},{\mathrm{C}}_{(2)},{\mathrm{C}}_{(3)},{\mathrm{C}}_{(4)},\mathrm{O},{\mathrm{N}}_{(2)},{\mathrm{N}}_{(3)}\}$; ${\mathrm{\Lambda }}_9=\{\mathrm{H},\mathrm{C},\mathrm{O},\mathrm{N},{\mathrm{S}}_{(2)},{\mathrm{S}}_{(4)},{\mathrm{S}}_{(6)},\mathrm{C}\text{l}\}$; ${\mathrm{\Lambda }}_{10}=\{\mathrm{H},{\mathrm{C}}_{(2)},{\mathrm{C}}_{(3)},{\mathrm{C}}_{(4)},{\mathrm{C}}_{(5)},\mathrm{O},{\mathrm{N}}_{(1)},{\mathrm{N}}_{(2)},{\mathrm{N}}_{(3)},\mathrm{F}\}$; and ${\mathrm{\Lambda }}_{11}=\{\mathrm{H},{\mathrm{C}}_{(2)},{\mathrm{C}}_{(3)},{\mathrm{C}}_{(4)},\mathrm{O},{\mathrm{N}}_{(2)},{\mathrm{N}}_{(3)},{\mathrm{S}}_{(2)},{\mathrm{S}}_{(4)},{\mathrm{S}}_{(6)},\mathrm{C}\text{l}\}$, where ${𝚎}_{(i)}$ for a chemical element $𝚎$ and an integer $i\ge 1$ means that a chemical element $𝚎$ with valence $i$.

- $|D_{\pi }|$: the size of data set $D_{\pi }$ over $\mathrm{\Lambda }$ for the property $\pi$.

- $\underset{¯}{n},\overline{n}$: the minimum and maximum values of the number $n(ℂ)$ of non-hydrogen atoms in compounds $ℂ$ in $D_{\pi }$.

- $\underset{¯}{a},\overline{a}$: the minimum and maximum values of $a(ℂ)$ for $\pi$ over compounds $ℂ$ in $D_{\pi }$.

- $|\mathrm{\Gamma }|$: the number of different edge-configurations of interior-edges over the compounds in $D_{\pi }$.

- $|ℱ|$: the number of non-isomorphic chemical rooted trees in the set of all 2-fringe-trees in the compounds in $D_{\pi }$.

- $K$: the number of descriptors in a feature vector $f(ℂ)$.

Stage 2. The newly defined feature function in our chemical model without suppressing hydrogen in Section 3 is used. We standardize the range of each descriptor and the range $\{t\in ℝ\mid \underset{¯}{a}\le t\le \overline{a}\}$ of property values $a(ℂ),ℂ\in D_{\pi }$.

Stage 3. For each chemical property $\pi$, we select a penalty value ${\lambda }_{\pi }$ in the Lasso function from 36 different values from 0 to 100 by conducting linear regression as a preliminary experiment.

We conducted an experiment in Stage 3 to evaluate the performance of the prediction function based on cross-validation. For a property $\pi$, an execution of a cross-validation consists of five trials of constructing a prediction function as follows. First partition the data set $D_{\pi }$ into five subsets $D^{(k)}$, $k\in [1,5]$ randomly. For each $k\in [1,5]$, the $k$-th trial constructs a prediction function ${\eta }^{(k)}$ by conducting a linear regression with the penalty term ${\lambda }_{\pi }$ using the set $D_{\pi }\backslash D^{(k)}$ as a training data set. We used scikit-learn version 0.23.2 with Python 3.8.5 for executing linear regression with Lasso function. For each property, we executed ten cross-validations and we show the median of test ${\mathrm{R}}^2({\eta }^{(k)},D^{(k)}),k\in [1,5]$ over all ten cross-validations. Recall that a subset of descriptors is selected in linear regression with Lasso function and let $K^{\prime }$ denote the average number of selected descriptors over all 50 trials. The running time per trial in a cross-validation was at most one second.

Table 1 shows the results on Stages 2 and 3, where we denote the following:

- ${\lambda }_{\pi }$: the penalty value in the Lasso function selected for a property $\pi$, where $a\mathrm{E}b$ means $a\times {10}^b$.

- $K^{\prime }$: the average of the number of descriptors selected in the linear regression over all 50 trials in ten cross-validations.

- test ${\mathrm{R}}^2$: the median of test ${\mathrm{R}}^2$ over all 50 trials in ten cross-validations.

Recall that the adjacency-configuration for leaf-edges was introduced as a new descriptor in this paper. Without including this new descriptor, the test ${\mathrm{R}}^2$ for property Vis was 0.790, that for Lumo was 0.799 and that for Mp was 0.796, while the test ${\mathrm{R}}^2$ for each of the other properties in Table 1 was almost the same.

From Table 1, we observe that a relatively large number of properties admit a good prediction function based on linear regression. The number $K^{\prime }$ of descriptors used in linear regression is considerably small for some properties. For example of property Vd, the four descriptors most frequently selected in the case of $\mathrm{\Lambda }=\{\mathrm{H},\mathrm{O},\mathrm{C},\mathrm{N}\}$ are the number of non-hydrogen atoms; the number of interior-vertices $v$ with ${\text{deg}}_{{\mathrm{C}}^{\text{int}}}(v)=1$; the number of fringe-trees r-isomorphic to the chemical rooted tree ${\psi }_1$ in Fig. 5; and the number of leaf-edges with adjacency-configuration $(\mathrm{O},\mathrm{C},2)$. The eight descriptors most frequently selected in the case of $\mathrm{\Lambda }=\{\mathrm{H},\mathrm{O},\mathrm{C},\mathrm{N},\mathrm{C}\text{l},{\mathrm{P}}_{(3)},{\mathrm{P}}_{(5)}\}$ are the number of non-hydrogen atoms; the number of interior-vertices $v$ with ${\text{deg}}_{{\mathrm{C}}^{\text{int}}}(v)=1$; the number of exterior-vertices $v$ with $\alpha (v)=\mathrm{C}\text{l}$; the number of interior-edges with edge-configuration ${\gamma }_i,i=1,2$, where ${\gamma }_1=(\mathrm{C2},\mathrm{C2},2)$ and ${\gamma }_2=(\mathrm{C3},\mathrm{C4},1)$; and the number of fringe-trees r-isomorphic to the chemical rooted tree ${\psi }_i,i=1,2,3$ in Fig. 5.

For the 18 properties listed in Table 1, we used ANN to construct prediction functions. For this purpose, we used our newly proposed feature vector and the experimental setup as explained in Tanaka et al. [25]. From these computation experiments, we observe that for the properties Hc, Vd, Ha, Hf, U0, Alpha and Cv, the test ${\mathrm{R}}^2$ scores of the prediction functions obtained by Lasso linear regression is at least 0.05 more than those obtained by ANN. For the properties OptR, Sl and SfT, the test ${\mathrm{R}}^2$ scores of the prediction functions obtained by ANN is at least 0.05 more than those obtained by Lasso linear regression. For the other properties, the test ${\mathrm{R}}^2$ scores obtained by Lasso linear regression and ANN are comparable.

Results on Phase 2. We used a set of seven instances $I_{\mathrm{a}}$, $I_{\mathrm{b}}^i,i\in [1,4]$, $I_{\mathrm{c}}$ and $I_{\mathrm{d}}$ based on seed graphs prepared by Shi et al. [26] to execute Stages 4 and 5 in Phase 2. We here present their seed graphs $G_{\mathrm{C}}$ (see Supplementary Material for the details of instances $I_{\mathrm{a}}$, $I_{\mathrm{b}}^i,i\in [1,4]$, $I_{\mathrm{c}}$ and $I_{\mathrm{d}}$). The seed graph $G_{\mathrm{C}}$ of instance $I_{\mathrm{a}}$ is illustrated in Fig. 4a. The seed graph $G_{\mathrm{C}}^1$ (resp., $G_{\mathrm{C}}^i,i=2,3,4$) of instances $I_{\mathrm{b}}^1$ and $I_{\mathrm{d}}$ (resp., $I_{\mathrm{b}}^i,i=2,3,4$) is illustrated in Fig. 6.

Instance $I_{\mathrm{c}}$ has been introduced by Shi et al. [26] in order to infer a chemical graph ${ℂ}^{†}$ such that the core of ${ℂ}^{†}$ is the same as the core of chemical graph ${ℂ}_A$: CID 24822711 in Fig. 7a and the frequency of each edge-configuration in the non-core of ${ℂ}^{†}$ is the same as that of chemical graph ${ℂ}_B$: CID 59170444 illustrated in Fig. 7b. This means that the seed graph $G_{\mathrm{C}}$ of $I_{\mathrm{c}}$ is the core of ${ℂ}_A$ which is indicated by a shaded area in Fig. 7a.

Instance $I_{\mathrm{d}}$ has been introduced by Shi et al. [26] in order to infer a monocyclic chemical graph ${ℂ}^{†}$ such that the frequency vector of edge-configurations in ${ℂ}^{†}$ is a vector obtained by merging those of two chemical graphs ${ℂ}_A$: CID 10076784 and ${ℂ}_B$: CID 44340250 illustrated in Fig. 7c,d, respectively.

Stage 4. We executed Stage 4 for five properties $\pi \in \{$Hc, Vd, OptR, IhcLiq, Vis$\mathrm{\}}$.

For the MILP formulation $ℳ(x,y;{𝒞}_1)$ in Section 4, we use the prediction function ${\eta }_{w,b}$ that attained the median test ${\mathrm{R}}^2$ in Table 1. We used CPLEX version 12.10 to solve an MILP in Stage 4. Tables 2,3,4,5,6 show the computational results of the experiment in Stage 4 for the five properties, where we denote the following:

- ${\underset{¯}{y}}^{*},{\overline{y}}^{*}$: lower and upper bounds ${\underset{¯}{y}}^{*}$, ${\overline{y}}^{*}\in ℝ$ on the value $a(ℂ)$ of a chemical graph $ℂ$ to be inferred;

- $\mathrm{\#}$v (resp., $\mathrm{\#}$c): the number of variables (resp., constraints) in the MILP in Stage 4;

- I-time: the time (sec.) to solve the MILP in Stage 4;

- $n$: the number $n\left({ℂ}^{†}\right)$ of non-hydrogen atoms in the chemical graph ${ℂ}^{†}$ inferred in Stage 4; and

- ${\mathrm{n}}^{\mathrm{int}}$: the number ${\mathrm{n}}^{\mathrm{int}}\left({ℂ}^{†}\right)$ of interior-vertices in the chemical graph ${ℂ}^{†}$ inferred in Stage 4;

- $\eta (f({ℂ}^{†}))$: the predicted property value $\eta (f({ℂ}^{†}))$ of the chemical graph ${ℂ}^{†}$ inferred in Stage 4.

From Tables 2,3,4,5,6 we observe that an instance with a large number of variables and constraints takes more running time than those with a smaller size in general. We solved all instances in this experiment with our MILP formulation in a few seconds to around 30 seconds.

Fig. 8a–e illustrate the chemical graphs ${ℂ}^{†}$ inferred from $I_{\mathrm{c}}$ with $({\underset{¯}{y}}^{*},{\overline{y}}^{*})=(13700,13800)$ of Hc, $I_{\mathrm{b}}^2$ with $({\underset{¯}{y}}^{*},{\overline{y}}^{*})=(21,22)$ of Vd, $I_{\mathrm{b}}^4$ with $({\underset{¯}{y}}^{*},{\overline{y}}^{*})=(70,71)$ of OptR, $I_{\mathrm{d}}$ with $({\underset{¯}{y}}^{*},{\overline{y}}^{*})=(1190,1210)$ of IhcLiq, and $I_{\mathrm{b}}^3$ with $({\underset{¯}{y}}^{*},{\overline{y}}^{*})=(1.85,1.90)$ of Vis, respectively.

Similarly, we executed Stage 4 for these seven instances $I_{\mathrm{a}}$, $I_{\mathrm{b}}^i,i\in [1,4]$, $I_{\mathrm{c}}$ and $I_{\mathrm{d}}$ for five properties $\pi \in \{$Hc, Vd, OptR, IhcLiq, Vis$\mathrm{\}}$ by using the prediction functions obtained by ANN. We list the running time to solve MILP formulation for each of these instances in Tables 7,8. From the computation experiments, we observe that for many instances, the running time is significantly faster than that of Stage 4 based on ANN.

Inferring a chemical graph with target values in multiple properties Once we obtained prediction functions ${\eta }_{\pi }$ for several properties $\pi$, include MILP formulations for these functions ${\eta }_{\pi }$ into a single MILP $ℳ(x,y;{𝒞}_1)$ so as to infer a chemical graph that satisfies given target values $y^{*}$ for these properties at the same time. As an additional experiment in Stage 4, we inferred a chemical graph that has a desired predicted value each of three properties Kow, Lp and Sl, where we used the prediction function ${\eta }_{\pi }$ for each property $\pi \in \{$Kow, Lp, Sl$\mathrm{\}}$ constructed in Stage 3. Table 9 shows the result of Stage 4 for inferring a chemical graph ${ℂ}^{†}$ from instances $I_{\mathrm{b}}^2$, $I_{\mathrm{b}}^3$ and $I_{\mathrm{b}}^4$ with $\mathrm{\Lambda }=\{\mathrm{H},\mathrm{C},\mathrm{N},\mathrm{O},{\mathrm{S}}_{(2)},{\mathrm{S}}_{(6)},\mathrm{Cl}\}$, where we denote the following:

- $\pi$: one of the three properties Kow, Lp and Slused in the experiment;

- ${\underset{¯}{y}}_{\pi }^{*},{\overline{y}}_{\pi }^{*}$: lower and upper bounds ${\underset{¯}{y}}_{\pi }^{*},{\overline{y}}_{\pi }^{*}\in ℝ$ on the predicted property value ${\eta }_{\pi }(f({ℂ}^{†}))$ of property $\pi \in \{$Kow, Lp, Sl$\mathrm{\}}$ for a chemical graph ${ℂ}^{†}$ to be inferred;

- $\mathrm{\#}$v (resp., $\mathrm{\#}$c): the number of variables (resp., constraints) in the MILP in Stage 4;

- I-time: the time (sec.) to solve the MILP in Stage 4;

- $n$: the number $n({ℂ}^{†})$ of non-hydrogen atoms in the chemical graph ${ℂ}^{†}$ inferred in Stage 4;

- $n^{\text{int}}$: the number $n^{\text{int}}({ℂ}^{†})$ of interior-vertices in the chemical graph ${ℂ}^{†}$ inferred in Stage 4; and

- ${\eta }_{\pi }(f({ℂ}^{†}))$: the predicted property value ${\eta }_{\pi }(f({ℂ}^{†}))$ of property $\pi \in \{$Kow, Lp, Sl$\mathrm{\}}$ for the chemical graph ${ℂ}^{†}$ inferred in Stage 4.

Fig. 8f illustrates the chemical graph ${ℂ}^{†}$ inferred from $I_{\mathrm{b}}^4$ with $({\underset{¯}{y}}_{{\pi }_1}^{*},{\overline{y}}_{{\pi }_1}^{*})=(-7.50,$$-7.40)$, $({\underset{¯}{y}}_{{\pi }_2}^{*},{\overline{y}}_{{\pi }_2}^{*})=(-0.70,-0.60)$ and $({\underset{¯}{y}}_{{\pi }_3}^{*},{\overline{y}}_{{\pi }_3}^{*})=(-11.4,-11.2)$ for ${\pi }_1=$Kow, ${\pi }_2=$Lp and ${\pi }_3=$Sl, respectively.

Stage 5. We executed Stage 5 to generate more target chemical graphs ${ℂ}^{*}$, where a chemical graph ${ℂ}^{*}$ is called a chemical isomer of a target chemical graph ${ℂ}^{†}$ of a topological specification $\sigma$ if $f({ℂ}^{*})=f({ℂ}^{†})$ and ${ℂ}^{*}$ also satisfies the same topological specification $\sigma$. We computed chemical isomers ${ℂ}^{*}$ of each target chemical graph ${ℂ}^{†}$ inferred in Stage 4. We executed an algorithm to generate chemical isomers of ${ℂ}^{†}$ up to 100 when the number of all chemical isomers exceeds 100. We can obtain such an algorithm from the dynamic programming proposed by Tanaka et al. [25] with a slight modification. The algorithm first decomposes ${ℂ}^{†}$ into a set of acyclic chemical graphs, next replaces each acyclic chemical graph $T$ with another acyclic chemical graph $T^{\prime }$ that admits the same feature vector as that of $T$, and finally assembles the resulting acyclic chemical graphs into a chemical isomer ${ℂ}^{*}$ of ${ℂ}^{†}$. Also, a lower bound on the total number of all chemical isomers of ${ℂ}^{†}$ can be computed by the algorithm without generating all of them.

Tables 2,3,4,5,6 show the computational results of the experiment in Stage 5 for the five properties, where we denote the following:

- D-time: the running time (sec.) to execute the dynamic programming algorithm in Stage 5 to compute a lower bound on the number of all chemical isomers ${ℂ}^{*}$ of ${ℂ}^{†}$ and generate all (or up to 100) chemical isomers ${ℂ}^{*}$;

- $ℂ$-LB: a lower bound on the number of all chemical isomers ${ℂ}^{*}$ of ${ℂ}^{†}$; and

- $\mathrm{\#}ℂ$: the number of all (or up to 100) chemical isomers ${ℂ}^{*}$ of ${ℂ}^{†}$ generated in Stage 5.

From Tables 2,3,4,5,6, we observe that for many cases the running time for generating up to 100 target chemical graphs in Stage 5 is less than 0.4 seconds. For some chemical graph ${ℂ}^{†}$, no chemical isomer was found by our algorithm. This is because each acyclic chemical graph in the decomposition of ${ℂ}^{†}$has no alternative acyclic chemical graph than the original one. On the other hand, some chemical graph ${ℂ}^{†}$ such as the one in $I_{\mathrm{d}}$ in Table 2 admits an extremely large number of chemical isomers ${ℂ}^{*}$. Remember that we know such a lower bound $ℂ$-LB on the number of chemical isomers without generating all of them.