# Query Structure Modeling for Inductive Logical Reasoning
Over Knowledge Graphs

Siyuan Wang1, Zhongyu Wei1,2, Meng Han3,
Zhihao Fan1, Haijun Shan4, Qi Zhang5, Xuanjing Huang5
1School of Data Science, Fudan University, China
2Research Institute of Intelligent and Complex Systems, Fudan University,
China
3Huawei Poisson Lab, China 4CEC GienTech Technology Co., Ltd, China
5School of Computer Science, Fudan University, China
<EMAIL_ADDRESS>
<EMAIL_ADDRESS><EMAIL_ADDRESS>
Corresponding author

###### Abstract

Logical reasoning over incomplete knowledge graphs to answer complex logical
queries is a challenging task. With the emergence of new entities and
relations in constantly evolving KGs, inductive logical reasoning over KGs has
become a crucial problem. However, previous PLMs-based methods struggle to
model the logical structures of complex queries, which limits their ability to
generalize within the same structure. In this paper, we propose a structure-
modeled textual encoding framework for inductive logical reasoning over KGs.
It encodes linearized query structures and entities using pre-trained language
models to find answers. For structure modeling of complex queries, we design
stepwise instructions that implicitly prompt PLMs on the execution order of
geometric operations in each query. We further separately model different
geometric operations (i.e., projection, intersection, and union) on the
representation space using a pre-trained encoder with additional attention and
maxout layers to enhance structured modeling. We conduct experiments on two
inductive logical reasoning datasets and three transductive datasets. The
results demonstrate the effectiveness of our method on logical reasoning over
KGs in both inductive and transductive settings. 111Codes are publicly
available at https://github.com/WangsyGit/InductiveLR.

##  1 Introduction

Figure 1:  Examples of inductive logical reasoning over KGs: testing queries
contain unseen entities and relations (in red) during training. Each query is
associated with an intrinsic logical structure and its natural language
interpretation.

Logical reasoning over knowledge graphs (KGs) aims to answer complex logical
queries given large-scale KGs Guu et al. (2015); Hamilton et al. (2018).
Recent years have witnessed increasing attention on logical reasoning over
widely used KGs such as Freebase Bollacker et al. (2008), Yago Suchanek et al.
(2007), NELL Carlson et al. (2010) and Wikidata Vrandečić and Krötzsch (2014).
With missing relations in the KG, it is challenging to deduce correct answers
to complex queries by traversing the graph. Previous work primarily focuses on
transductive logical reasoning where the training and testing are done on the
same KG with the same group of entities. They typically rely on geometric
embedding-based methods to map both entities and queries into a joint low-
dimensional vector space Hamilton et al. (2018); Ren et al. (2020); Ren and
Leskovec (2020). The goal is to push the embeddings of answer entities and
queries to be close to each other, allowing answers to be predicted through
embedding similarity even when the involved relation is absent. In contrast,
the inductive setting of logical reasoning has been rarely studied which
requires generalizing to unseen entities and relations or even new KGs. As
real-world KGs are usually dynamic with emerging unseen entities and
relations, it’s significant to explore the inductive setting for complex query
answering.

Figure 2:  The overall architecture of our method for inductive KG logical
reasoning with structure knowledge modeling. [inter] and [proj] are
respectively abbreviations of [intersection] and [projection].

Existing research on inductive logical reasoning mainly follows two
directions. The first inherits embedding-based methods and incorporates type
as additional information to improve inductive capability Hu et al. (2022),
which can not generalize to unseen types of entities and relations. The second
direction leverages pre-trained language models (PLMs) to encode textual
information of entities/relations for generalization to unseen elements Wang
et al. (2021b); Daza et al. (2021); Wang et al. (2021a). PLMs-based approaches
provide more flexible solutions and generate better results. However, they
only explore link prediction tasks of one-step reasoning, and simply linearize
the triplet or subgraph into text sequence without modeling explicit reasoning
structure Yao et al. (2019); Zha et al. (2022). An example is shown in Figure
1. Two findings stand out. (1) The query $q_{1}$ and $q_{2}$ appear to be
similar in format (both as a conjunction of three terms) but actually have
different logical structures. PLMs-based methods that encode flattened queries
can not model this structure information for correct logical reasoning. (2)
Although queries $q_{1}$ and $q_{3}$ (also $q_{2}$ and $q_{4}$) contain
different elements, they share the same logical structure. Motivated by these,
we argue that structure modeling of different complex queries can further
boost the generalization ability of logical reasoners.

In this paper, we propose to model query structure for inductive logical
reasoning over KGs. Specifically, we transform the query structure into a
sequence using textual names of involved entities, relations, and logical
operators. For complex query structures composed of multiple geometric
operations over entities and relations, we introduce two measures to enable
logical structure modeling during text encoding. First, we design stepwise
instructions for different query types to indicate which operation in the
query structure should be conducted at each step and feed them as the
structural prompt to PLMs. Besides, we extend the pre-trained encoder with an
additional attention layer and a maxout layer to respectively model different
geometric operations including projection, intersection, and union on the
representation space, to implicitly inject structured modeling into PLMs. Our
proposed method is a generic inductive framework, which can be plugged into
different PLMs for better performance.

We conduct experiments on two datasets for inductive logical reasoning over
KGs, FB15k-237-V2 and NELL-V3 Teru et al. (2020) as well as three transductive
datasets, FB15k Bordes et al. (2013), FB15k-237 Toutanova and Chen (2015), and
NELL995 Xiong et al. (2017). The results demonstrate that our method achieves
strong inductive performance on unseen entities and relations, even across
different KGs, without sacrificing logical reasoning capability and
generalizability to new query structures.

##  2 Methodology

In this work, we study the task of complex logical reasoning over KGs. The
input is a first-order logic query $q$ which can include any set of
existential quantification ($\exists$), conjunction ($\wedge$), and
disjunction ($\vee$) operators (such as the query in Figure 1). Our goal is to
predict a set of entities $\mathcal{A}=\\{a_{1},a_{2},a_{3}...\\}$ that answer
the query $q$ based on an incomplete KG $\mathcal{G=(E,R)}$ which consists of
a set of triplets $(h,r,t)$ but lacks several involved relations. Here
$h,t\in\mathcal{E}$ are the head and tail entities and $r\in\mathcal{R}$ is
the relation between them. We mainly focus on the inductive setting of KG
logical reasoning, where the evaluated queries contain entities/relations that
are completely unseen during the training period.

Figure 2 shows the overall architecture of our model. We propose to encode the
text sequences of query structures and predict the answer entities based on
representation similarity for inductive logical reasoning over KGs. In this
section, we first list different types of query structures studied in logical
reasoning over KGs(§ 2.1). Then according to various query types, we introduce
our structure linearization and structural prompts for textual encoding(§
2.2). The geometric operation modeling and query answer prediction modules are
described in (§ 2.3) and (§ 2.4). Finally, we provide the details about
training and inference (§ 2.5).

###  2.1 Query Structure Types

Following Ren et al. (2020), we consider 9 types of complex query structures
which are composed of different sets of geometric operations (including
projection, intersection and union) over entities and relations. These include
six single-operation query structures and three mixed-operation ones.
Specifically, three query types only focus on projection, including one-
relation projection (1p), two-relation projection (2p), and three-relation
projection (3p). Two query types focus on the intersection of two triplets
(2i) and three triplets (3i), and another one focuses on the union of two
triplets (2u). The three mixed-operation query structures are respectively the
combinations of intersection&projection (ip), projection&intersection (pi),
and union&projection (up). The different query structures are illustrated as
the following formula:

$\displaystyle[1p]\ q=v_{?}$ $\displaystyle:r_{1}(e_{1},v_{?})$
$\displaystyle[2p]\ q=v_{?}$ $\displaystyle.\exists v:r_{1}(e_{1},v)\wedge
r_{2}(v,v_{?})$ $\displaystyle[3p]\ q=v_{?}$ $\displaystyle.\exists
v_{1},v_{2}:r_{1}(e_{1},v_{1})\wedge r_{2}(v_{1},v_{2})$ $\displaystyle\wedge
r_{3}(v_{2},v_{?})$ $\displaystyle[2i]\ q=v_{?}$
$\displaystyle:r_{1}(e_{1},v_{?})\wedge r_{2}(e_{2},v_{?})$
$\displaystyle[3i]\ q=v_{?}$ $\displaystyle:r_{1}(e_{1},v_{?})\wedge
r_{2}(e_{2},v_{?})\wedge r_{3}(e_{3},v_{?})$ $\displaystyle[pi]\ q=v_{?}$
$\displaystyle.\exists v:r_{1}(e_{1},v)\wedge r_{2}(v,v_{?})\wedge
r_{3}(e_{2},v_{?})$ $\displaystyle[ip]\ q=v_{?}$ $\displaystyle.\exists
v:r_{1}(e_{1},v)\wedge r_{2}(e_{2},v)\wedge r_{3}(v,v_{?})$
$\displaystyle[2u]\ q=v_{?}$ $\displaystyle:r_{1}(e_{1},v_{?})\vee
r_{2}(e_{2},v_{?})$ $\displaystyle[up]\ q=v_{?}$ $\displaystyle.\exists
v:(r_{1}(e_{1},v)\vee r_{2}(e_{2},v))\wedge r_{3}(v,v_{?})$ (1)

where $e_{i}$ and $v_{i}$ are the anchor entities and existentially quantified
bound variables entities, and $v_{?}$ are the target answer entities to the
query. As these complex queries contain rich logical structural information,
we need to model the structure knowledge during textual encoding for better
inductive generalization within the same logical structure.

###  2.2 Query Structure Encoding

In order to use PLMs for better generalization to unseen entities/relations,
we first need to linearize the query structures into text sequences. We also
design instructions for each query type as a structural prompt to implicitly
indicate the order of geometric operations execution to PLMs. We concatenate
the linearized query and structural prompt as the input, and encode them to
obtain the query representation for matching with the answer entities.

#### Query Structure Linearization

Given a query structure $q$, we customize its linearization method according
to the query type. For each triplet $r_{i}(e_{i},v)$ in the query, we
formulate it as “[anchor] $t(e_{i})$ [relation] $t(r_{i})$” and ignores the
intermediate variable entities, where $t(e_{i})$ and $t(r_{i})$ are the
textual names of anchor entity $e_{i}$ and relation $r_{i}$. Then we add the
names of logical operations before the involved subqueries. For example, the
query structure of type [2p] can be linearized into sequence “[projection]
[anchor] $t(e_{1})$ [relation] $t(r_{1})$ [projection] [relation] $t(r_{2})$”
and the query structure of type [2i] can be mapped into “[intersection]
[projection] [anchor] $t(e_{1})$ [relation] $t(r_{1})$ [projection] [anchor]
$t(e_{2})$ [relation] $t(r_{2})$”.

For query types [ip] and [up] that are composed of intersection/union and
projection, the last relation projection is conducted over the
intersection/union of previous triplets. Directly flattening the query is
unable to keep such structural information. Therefore, we propose to split the
intersection/union and repeatedly connect each triplet with the last relation
projection which moves the intersection/union operation to the final step.
This transformation is equivalent to the original query structure. For
example, the following two query structures are equivalent and are both of
type [up].

$\displaystyle(r_{1}(e_{1},v)\vee r_{2}(e_{2},v))\wedge r_{3}(v,v_{?})$
$\displaystyle(r_{1}(e_{1},v)\wedge r_{3}(v,v_{?}))\vee(r_{2}(e_{2},v)\wedge
r_{3}(v,v_{?}))$ (2)

Based on this transformation, we linearize the query structure of type [up]
into the text sequence as “[union] [projection] [anchor] $t(e_{1})$ [relation]
$t(r_{1})$ [projection] [relation] $t(r_{3})$ [projection] [anchor] $t(e_{2})$
[relation] $t(r_{2})$ [projection] [relation] $t(r_{3})$”. The details of
structure linearization templates for each query type are listed in Appendix
A.

#### Structural Prompt

Besides feeding the linearized query structure into PLMs, we also introduce
stepwise instructions that indicate the order of geometric operations
execution to prompt the pre-trained encoder with implicit structural
information of the query. Specifically, each prompt consists of two parts: the
number of total execution steps and the operation order in the query, which is
formulated as “total steps: operation order”. For query types [ip] and [up],
we define the total steps and operation order according to the transformed
query structure in Eq 2.2. The detailed structural prompts of each query type
are presented in Table 1.

Query Type |  Structural Prompt
---|---
[1p] |  One step: [proj]
[2p] |  Two steps: [proj], then [proj]
[3p] |  Three steps: [proj], then [proj], then [proj]
[2i] |  Three steps: [proj], and [proj], [inter]
[3i] |  Four steps: [proj], and [proj], and [proj], [inter]
[pi] |  Four steps: [proj], then [proj], and [proj], [inter]
[ip] |  Five steps: [proj], then [proj], and [proj], then [proj], [inter]
[2u] |  Three steps: [proj], and [proj], [union]
[up] |  Five steps: [proj], then [proj], and [proj], then [proj], [union]

Table 1:  Structural prompts of sifferent query types. “[proj]” and “[inter]”
are the abbreviations of “[projection]” and “[intersection]”.

Then the input $t$ can be obtained by concatenating the structural prompt $s$
and the linearized query structure $t(q)$, which is formulated as “[CLS]
[qtype] s [SEP] t(q)”. We feed it to the pre-trained encoder and obtain the
output hidden states $H=(h_{1},h_{2},...,h_{|t|})$.

###  2.3 Geometric Operation Modeling

To further enhance the structure modeling during textual encoding, we propose
to separately model different geometric operations in logical queries to
explore their spatial characteristics. As Transformer-based encoders are
widely used to implicitly learn the translation function for simple link
prediction and question answering Yao et al. (2019); Wang et al. (2021a), we
directly utilize it for modeling multi-hop relation projection. For each path
of relation projection $r_{1}(e_{1},v_{1})\wedge...\wedge r_{j}(v_{j},v_{?})$,
we extract the hidden states corresponding to sequence “[anchor] $t(e_{1})$
[relation] $t(r_{1})$ ... [projection] [relation] $t(r_{j})$” from
$H=(h_{1},h_{2},...,h_{|t|})$. We then take the average as the representation
of target entity $v_{?}$, which also can be viewed as the representation
$h_{q}$ of the query that only involves relation projection operation.

For the intersection and union of multiple subqueries, we adopt an attention
layer Bahdanau et al. (2014) and a maxout layer Goodfellow et al. (2013) on
top of the pre-trained encoder to respectively model these two operations in
the geometric representation space. Specifically, we feed the representations
of target entities in all subqueries to these two additional layers to achieve
the intersection and union operations. The output can be taken as the query
representation $h_{q}$ that contains intersection or union.

As presented in Ren and Leskovec (2020), the complete set of first-order logic
operations encompasses existential quantification ($\exists$), conjunction
($\wedge$), disjunction ($\vee$) and negation($\neg$). Our approach covers the
first three operations by modeling relation projection, intersection, and
union respectively. The negation operation is not individually modeled as pre-
trained encoders are capable of capturing semantic exclusion for negation. We
can add negative terms such as “not” before the corresponding relations within
the input and feed it into pre-trained encoders to naturally address this
task.

###  2.4 Query Answer Prediction

To answer the complex query $q$, we adopt the Siamese dual encoder Gillick et
al. (2018) to respectively encode the query $q$ and the candidate entity
$c_{i}$ to match the answer. We formulate the entity input as the sequence
“[CLS] [target] $t(c_{i})$” and feed it into the pre-trained encoder to obtain
the candidate entity representation $h_{c_{i}}$ by taking the average of the
hidden states. Then we compute the similarity $d_{i}$ between the query
representation $h_{q}$ and entity representation $h_{c_{i}}$, and encourage
the query representation to be similar to the positive answer entities while
dissimilar to negative entities. The entities whose representations are
similar enough to the query will be predicted as the answers. We can pre-
compute the representations of all candidate entities and utilize them to
predict answers for different queries more efficiently.

The above matching scheme can handle the inductive setting when the candidate
entity set is not closed and new entities may arise. To improve the answer
prediction accuracy in the transductive setting where the candidate entity set
is closed, we also employ a classification layer on top of the query
representation $h_{q}$. Given the fixed candidate entity set
$(c_{1},c_{2},...,c_{N})$, the classification layer with a softmax function
will output an $N$-dimensional plausibility distribution
$(s_{c_{1}},s_{c_{2}},...,s_{c_{N}})$ for each candidate entity $c_{i}$ to be
the answer.

###  2.5 Training & Inference

We simultaneously optimize a matching objective and a classification objective
to train our inductive model for answering complex logical queries. For the
former, we adopt contrastive learning Chen et al. (2020) which needs to
separate the positive and negative answers for the query. We take the given
ground-truth answer $c_{+}$ as positive and implement in-batch negative
sampling to collect the negatives. We measure the similarity between the query
and entities using dot product, and follow He et al. (2020) to utilize InfoNCE
as the contrastive loss. The loss function is formulated as Eq. 3 where $\tau$
is

$\displaystyle\mathcal{L}_{M}$ $\displaystyle=-\log\frac{\exp(h_{q}\cdot
h_{c_{+}}/\tau)}{\sum_{i=1}^{N}(\exp(h_{q}\cdot h_{c_{i}}/\tau)}$ (3)

the temperature hyper-parameter and $N$ is the total number of candidate
entities including positives and negatives. For the classification objective,
we take all entities in each KG as candidate entities and calculate the cross-
entropy loss as Eq. 4.

$\displaystyle\mathcal{L}_{C}=-\log\frac{s_{c_{+}}}{\sum_{i=1}^{N}\exp(s_{c_{i}})}$
(4)

These two losses are combined in a weighted manner as
$\mathcal{L}=\mathcal{L}_{M}+\lambda\mathcal{L}_{C}$ and $\lambda$ is the
weighted hyper-parameter.

During inference, we perform differently for inductive and transductive
logical query answering. For the inductive reasoning, we utilize the matching
scheme and rank the representation similarities between the query and all
candidate entities for query answer prediction. For the transductive
inference, we only adopt the classification scheme and find the most plausible
answer according to classification scores.

##  3 Experiments

###  3.1 Experiment Setup

We conduct experiments on complex logical reasoning over KGs in both inductive
and transductive setting. For the inductive setting, we adopt two datasets,
FB15k-237-V2 and NELL-V3, that have disjoint sets of entities for training and
evaluation, as introduced by Teru et al. (2020). To further challenge our
model, we also illustrate the cross-KG inductive generalization performance by
respectively taking FB15k and NELL995 for training/inference and
inference/training that contain completely different entities and relations.
In the transductive setting, we evaluate our model on the generated queries
from three datasets: FB15k Bordes et al. (2013), FB15k-237 Toutanova and Chen
(2015), and NELL995 Xiong et al. (2017), as proposed by Ren et al. (2020). All
these datasets cover nine types of query structures. We follow the setting of
Ren et al. (2020) to illustrate the generalization within the same structure
to unseen entities and relations, and also the generalization to more
complicated unseen structures composed of different structures. Specifically,
we train our model on the first five types of query structures (1p, 2p, 3p,
2i, 3i) and evaluate it on all nine query types (1p, 2p, 3p, 2i, 3i, pi, ip,
2u, up), including both seen and unseen query structures during training.

Model | Avg | 1p | 2p | 3p | 2i | 3i | ip | pi | 2u | up
---|---|---|---|---|---|---|---|---|---|---
FB15k-237-V2
_Q2B_ | 0.043 | 0.005 | 0.055 | 0.017 | 0.007 | 0.007 | 0.078 | 0.129 | 0.027 | 0.061
_TEMP(GQE)_ | 0.163 | 0.146 | 0.221 | 0.141 | 0.139 | 0.144 | 0.157 | 0.220 | 0.097 | 0.201
_BiQE_ | 0.158 | 0.286 | 0.151 | 0.107 | 0.187 | 0.240 | 0.125 | 0.155 | 0.085 | 0.091
_SILR_ | 0.178 | 0.309 | 0.121 | 0.106 | 0.237 | 0.274 | 0.148 | 0.181 | 0.110 | 0.116
NELL-V3
_Q2B_ | 0.017 | 0.002 | 0.022 | 0.005 | 0.003 | 0.002 | 0.026 | 0.048 | 0.018 | 0.028
_TEMP(GQE)_ | 0.062 | 0.096 | 0.057 | 0.060 | 0.072 | 0.081 | 0.047 | 0.062 | 0.040 | 0.040
_BiQE_ | 0.089 | 0.178 | 0.081 | 0.081 | 0.082 | 0.092 | 0.067 | 0.081 | 0.066 | 0.069
_SILR_ | 0.101 | 0.197 | 0.079 | 0.074 | 0.103 | 0.122 | 0.094 | 0.090 | 0.080 | 0.068

Table 2:  Inductive H@10 results of different structured queries on
FB15k-237-V2 and NELL-V3 datasets. The results of _Q2B_ and _TEMP(Q2B)_ are
taken from Hu et al. (2022). Avg is the average performance of all query
types.

Model | Avg | 1p | 2p | 3p | 2i | 3i | ip | pi | 2u | up
---|---|---|---|---|---|---|---|---|---|---
FB15k $\rightarrow$ NELL995
_BiQE_ | 0.042 | 0.077 | 0.076 | 0.054 | 0.025 | 0.023 | 0.036 | 0.035 | 0.016 | 0.039
_SILR_ | 0.069 | 0.166 | 0.101 | 0.096 | 0.033 | 0.040 | 0.020 | 0.027 | 0.063 | 0.073
NELL995 $\rightarrow$ FB15k
_BiQE_ | 0.082 | 0.175 | 0.103 | 0.078 | 0.082 | 0.094 | 0.041 | 0.057 | 0.033 | 0.072
_SILR_ | 0.098 | 0.218 | 0.128 | 0.100 | 0.081 | 0.089 | 0.029 | 0.059 | 0.115 | 0.065

Table 3:  Inductive H@10 results on cross-KG generalization. FB15k
$\rightarrow$ NELL995 and NELL995 $\rightarrow$ FB15k respectively mean that
the model is trained on FB15k and tested on NELL995, and vice versa.

The data split statistics of logical queries in these datasets are provided in
Appendix B.

###  3.2 Implementation Details

We take bert-large-cased and bert-base-cased Devlin et al. (2018) as the pre-
trained encoder for encoding the query structure in (FB15k-237-V2, NELL-V3)
and (FB15k, FB15k-237, and NELL995), respectively. All models are implemented
using Huggingface Wolf et al. (2019), and trained for 30 epochs on 4 NVIDIA
Tesla V100 GPUs with 16 GB of memory. The Adam is taken as the optimizer and
the learning rate is 1.5e-4. We use a linear learning rate scheduler with
$10\%$ warmup proportion. The weight hyper-parameter to balance losses is set
to $\lambda=0.3$ or $\lambda=0.4$. For automatic evaluation, we use Hits@$K$
($K=3,10$) as the metrics, which calculate the proportion of correct answer
entities ranked among the top-$K$.

Model | Avg | 1p | 2p | 3p | 2i | 3i | ip | pi | 2u | up
---|---|---|---|---|---|---|---|---|---|---
FB15k-237
_GQE_ | 0.228 | 0.402 | 0.213 | 0.155 | 0.292 | 0.406 | 0.083 | 0.17 | 0.169 | 0.163
_Q2B_ | 0.268 | 0.467 | 0.240 | 0.186 | 0.324 | 0.453 | 0.108 | 0.205 | 0.239 | 0.193
_TEMP(Q2B)_ | 0.294 | 0.457 | 0.278 | 0.234 | 0.369 | 0.496 | 0.229 | 0.117 | 0.276 | 0.189
_BiQE_ | 0.262 | 0.472 | 0.293 | 0.245 | 0.345 | 0.473 | 0.096 | 0.211 | 0.075 | 0.145
_SILR_ | 0.296 | 0.471 | 0.302 | 0.249 | 0.358 | 0.484 | 0.113 | 0.222 | 0.283 | 0.181
FB15k
_GQE_ | 0.386 | 0.636 | 0.345 | 0.248 | 0.515 | 0.624 | 0.151 | 0.310 | 0.376 | 0.273
_Q2B_ | 0.484 | 0.786 | 0.413 | 0.303 | 0.593 | 0.712 | 0.211 | 0.397 | 0.608 | 0.330
_Q2B+TEMP_ | 0.554 | 0.840 | 0.498 | 0.422 | 0.674 | 0.779 | 0.483 | 0.261 | 0.727 | 0.300
_BiQE_ | 0.482 | 0.810 | 0.526 | 0.461 | 0.677 | 0.781 | 0.209 | 0.456 | 0.185 | 0.232
_SILR_ | 0.596 | 0.867 | 0.575 | 0.489 | 0.736 | 0.824 | 0.287 | 0.500 | 0.767 | 0.321
NELL995
_GQE_ | 0.247 | 0.418 | 0.228 | 0.205 | 0.316 | 0.447 | 0.081 | 0.186 | 0.199 | 0.139
_Q2B_ | 0.306 | 0.555 | 0.266 | 0.233 | 0.343 | 0.480 | 0.132 | 0.212 | 0.369 | 0.163
_Q2B+TEMP_ | 0.373 | 0.625 | 0.343 | 0.342 | 0.410 | 0.552 | 0.209 | 0.141 | 0.477 | 0.262
_BiQE_ | 0.309 | 0.632 | 0.310 | 0.332 | 0.370 | 0.525 | 0.091 | 0.164 | 0.173 | 0.184
_SILR_ | 0.342 | 0.641 | 0.329 | 0.337 | 0.376 | 0.532 | 0.055 | 0.177 | 0.450 | 0.182

Table 4:  Transductive H@3 results on FB15k-237, FB15k and NELL995 datasets.
The results of _GQE_ and _Q2B_ are taken from Ren et al. (2020).

###  3.3 Inductive Setting

#### Unseen Entities Generalization

To illustrate the inductive performance of complex logical reasoning over KGs,
we first make a comparison on FB15k-237-V2 and NELL-V3 datasets for
generalizing to unseen entities. We compare our structure-modeled inductive
logical reasoning method (_SILR_), with the baseline embedding-based method
_Q2B_ Ren et al. (2020) and the best version of the inductive model
_TEMP(GQE)_ Hu et al. (2022). _TEMP(GQE)_ enriches embedding method _GQE_
Hamilton et al. (2018) with type information which has achieved the state-of-
the-art performance. _BiQE_ Kotnis et al. (2021) is also a textual encoding
method with positional embedding for logical query but only in the
transductive setting. We reimplement it by replacing the original
classification-based prediction with the matching scheme for an inductive
comparison.

The experimental results are shown in Table 2. We can see that our _SILR_
outperforms all other models on both FB15k-237-V2 and NELL-V3 datasets by a
considerable margin. This highlights the effectiveness of our method for
inductive logical reasoning over unseen entities. Additionally, the
improvement over the other positional textual encoding model, _BiQE_ ,
demonstrates that our structure knowledge modeling during textual encoding is
capable of enhancing the inductive complex query answering capability.

#### Cross-KG Generalization

We further explore a more challenging cross-KG inductive setting, where the
model is trained and tested on different datasets and requires generalizing to
completely different KGs. Specifically, we take FB15k and NELL995 as the
source/target and target/source datasets, respectively. In this scenario, we
adopt the few-shot setting, where 500 random samples of the target domain are
provided for continual learning to achieve better transferring. As embedding-
based methods, even when aware of type information, are unable to embed most
entities and relations in new KGs with unseen types, we only compare our
_SILR_ with the reimplemented _BiQE_. The results in Table 3 show that our
_SILR_ performs better than _BiQE_ , and it can not only generalize to unseen
entities but also perform logical reasoning over new KGs with only a few
portions observed. This manifests the effectiveness of our method on inductive
logical reasoning over KGs even in the real-world challenging cross-KG
setting.

###  3.4 Transductive Setting

Although textual encoding methods have the inductive potential, their
performance often lags behind embedding-based models due to learning
inefficiency and the inability to structure knowledge modeling Wang et al.
(2022). We also compare our _SILR_ with transductive logical reasoning methods
to illustrate the logical reasoning performance over KGs with structure
modeling. The compared models including _GQE_ , _Q2B_ , _TEMP(Q2B)_ and _BiQE_
where the first three are embedding-based models while the last one is a
positional textual encoding model. Since _BiQE_ does not evaluate query types
2u and up, and does not provide results for FB15k, we reimplement it for a
fair comparison.

The results are shown in Table 4. Our _SILR_ outperforms _BiQE_ , particularly
in query types involving intersection and union, indicating that our structure
knowledge modeling can effectively improve the logical reasoning performance
of textual encoding and help generalize to unseen query structures. Although
textual encoding methods have the potential for inductive KG reasoning, they
still lag behind embedding-based methods for the transductive setting, due to
explicit structure modeling and better learning efficiency of embedding-based
methods Wang et al. (2022). In this work, we mainly focus on improving textual
encoding methods for inductive complex reasoning, but our method still
achieves comparable transductive performance. This demonstrates the
effectiveness of our inductive method with query structure modeling on
transductive logical reasoning over KGs.

###  3.5 Further Analysis

#### Ablation Study

To dive into the impact of different components in our model on both inductive
and transductive logical reasoning over KGs, we conduct an ablation study on
the FB15k-237-V2 and FB15k-237 datasets. We respectively take bert-large-cased
and bert-base-cased as baseline models for FB15k-237-V2 and FB15k-237. They
remove Structural Prompt (_SP_) and Geometric Operation Modeling (_GOM_) from
the final model _SILR_ , which directly encodes linearized query structures
for answer prediction.

As shown in Table 5, incorporating structural prompting and geometric
operation modeling can both improve the baseline model, but still perform
worse than our final _SILR_. This indicates that these two measures for
modeling structure knowledge during query text encoding can enhance the
inductive and transductive performance of logical reasoning over KGs.

Model | FB15k-237-V2 | FB15k-237
---|---|---
H@3 | H@10 | H@3 | H@10
_Baseline_ | 0.072 | 0.151 | 0.258 | 0.397
 _w/ SP_ | 0.073 | 0.162 | 0.288 | 0.429
 _w/ GOM_ | 0.071 | 0.167 | 0.281 | 0.420
_SILR_ | 0.079 | 0.178 | 0.296 | 0.437

Table 5:  Ablation study results on FB15k-237-V2 and FB15k-237. _SP_ and _GOM_
respectively denote Structural Prompt and Geometric Operation Modeling.

#### Query Structure Generalization

Embedding-based methods are known to generalize well to unseen query
structures due to their explicit spatial structure learning. To analyze the
generalizability of our implicit structure-modeled textual encoding method to
different logical query structures, we further construct four types of
complicated query structures with more relations and geometric operations,
including 4p, 5p, 3ip and i2p, based on the query structures in test sets. The
detailed illustrations and explanations of these query structures are given in
Appendix C. We directly evaluate our method on these more complicated queries
in both inductive and transductive datasets FB15k-237-V2 and FB15k-237. The
results are listed in Table 6. We can see that compared to seen and unseen
query structures in the original datasets in Table 2 and 4, our method can
also generalize to these complicated structures with more logical operations
and achieves impressive performance. This demonstrates that the design of
structural prompt and implicit geometric operation modeling can effectively
learn structure knowledge and improve query structure generalizability.

Model | 4p | 5p | 3ip | i2p
---|---|---|---|---
Inductive | 0.186 | 0.089 | 0.256 | 0.092
Transductive | 0.192 | 0.188 | 0.367 | 0.199

Table 6:  Analysis of generalization to more complicated query structures. For
inductive, we report the H@10 scores on FB15k-237-V2 and for transductive we
report the H@3 results on FB15k-237.

##  4 Related Work

Answering first-order logic queries over incomplete KGs is a challenging task
Guu et al. (2015). Previous research Lin et al. (2015); Hamilton et al. (2018)
mainly studies transductive logical reasoning, where the training and testing
are performed on the same KG. Embedding-based methods are widely used to embed
both logical queries and entities into a joint low-dimensional vector space
and push answer entities and queries to be close enough, enabling answer
prediction through embedding similarity, even when the involved relation is
absent. Following this paradigm, some further propose extending the embedding
of the query/entity from a single point to a region Ren et al. (2020); Zhang
et al. (2021); Choudhary et al. (2021b) or probabilistic distribution Ren and
Leskovec (2020); Choudhary et al. (2021a) over vector space to map arbitrary
first-order logic queries into sets of answer entities. However, these methods
are unable to tackle the inductive problem, which requires generalizing to
unseen entities and relations. Although Hu et al. (2022) proposes enriching
the entity and relation embedding with type information for inductive logical
reasoning, it can only generalize to elements of observed types.

Another line of research focuses on inductive logical reasoning over KGs using
textual encoding methods. With the advance of large-scale pre-trained language
models Devlin et al. (2018); Liu et al. (2019), these methods propose
transforming the graph structures into linearized text and utilize PLMs for
encoding Yao et al. (2019); Zha et al. (2022). With the strong
generalizability of PLMs, they can easily generalize to unseen
entities/relations, but struggle to model structure knowledge during text
encoding. Some works Wang et al. (2021b); Daza et al. (2021); Wang et al.
(2021a) propose to follow TransE Bordes et al. (2013) to apply the translation
function between entity and relation representations for geometric structure
learning. Nevertheless, these methods usually require descriptions of entities
and relations for encoding and assume these descriptions are readily
available. Besides, they only focus on simple link prediction tasks without
exploring complex structure modeling in logical reasoning, which is essential
for generalizing within the same query structure/type. We thus propose to
simultaneously encode the linearized query and preserve the logical structure
knowledge by structural prompt and separate geometric operation modeling for
inductive logical reasoning over KGs.

##  5 Conclusion

In this paper, we present the first flexible inductive method for answering
complex logical queries over incomplete KGs which can generalize to any unseen
entities and relations. To accomplish this goal, we propose a structure-model
textual encoding model that utilizes PLMs to encode linearized query
structures to find the answer entities. For structure modeling of complex
queries, we design structural prompts to implicitly indicate PLMs the order of
geometric operations execution in each query, and separately model three
geometric operations on representation space using a pre-trained encoder, an
attention layer, and a maxout layer. Experimental results demonstrate the
effectiveness of our model on logical reasoning over KGs in both inductive and
transductive settings.

## Limitations

This study has potential limitations. First, it only focuses on answering
existential positive first-order logic queries but does not support the
negation operation. We will later address this limitation by modeling the
negation operation. Second, we utilize BERT as the backbone model for
inductive generalization due to computing resource limits. We plan to
investigate the use of more powerful pre-trained language models with stronger
generalizability in future research to improve inductive logical reasoning
over KGs.

## Acknowledgments

This work is supported by National Natural Science Foundation of China (No.
6217020551) and Science and Technology Commission of Shanghai Municipality
Grant (No.21QA1400600).

## References

  * Bahdanau et al. (2014) Dzmitry Bahdanau, Kyunghyun Cho, and Yoshua Bengio. 2014.  Neural machine translation by jointly learning to align and translate.  _arXiv preprint arXiv:1409.0473_.
  * Bollacker et al. (2008) Kurt Bollacker, Colin Evans, Praveen Paritosh, Tim Sturge, and Jamie Taylor. 2008\.  Freebase: a collaboratively created graph database for structuring human knowledge.  In _Proceedings of the 2008 ACM SIGMOD international conference on Management of data_ , pages 1247–1250.
  * Bordes et al. (2013) Antoine Bordes, Nicolas Usunier, Alberto Garcia-Duran, Jason Weston, and Oksana Yakhnenko. 2013.  Translating embeddings for modeling multi-relational data.  _Advances in neural information processing systems_ , 26.
  * Carlson et al. (2010) Andrew Carlson, Justin Betteridge, Bryan Kisiel, Burr Settles, Estevam R Hruschka, and Tom M Mitchell. 2010.  Toward an architecture for never-ending language learning.  In _Twenty-Fourth AAAI conference on artificial intelligence_.
  * Chen et al. (2020) Ting Chen, Simon Kornblith, Mohammad Norouzi, and Geoffrey Hinton. 2020.  A simple framework for contrastive learning of visual representations.  In _International conference on machine learning_ , pages 1597–1607. PMLR.
  * Choudhary et al. (2021a) Nurendra Choudhary, Nikhil Rao, Sumeet Katariya, Karthik Subbian, and Chandan Reddy. 2021a.  Probabilistic entity representation model for reasoning over knowledge graphs.  _Advances in Neural Information Processing Systems_ , 34:23440–23451.
  * Choudhary et al. (2021b) Nurendra Choudhary, Nikhil Rao, Sumeet Katariya, Karthik Subbian, and Chandan K Reddy. 2021b.  Self-supervised hyperboloid representations from logical queries over knowledge graphs.  In _Proceedings of the Web Conference 2021_ , pages 1373–1384.
  * Daza et al. (2021) Daniel Daza, Michael Cochez, and Paul Groth. 2021.  Inductive entity representations from text via link prediction.  In _Proceedings of the Web Conference 2021_ , pages 798–808.
  * Devlin et al. (2018) Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova. 2018.  Bert: Pre-training of deep bidirectional transformers for language understanding.  _arXiv preprint arXiv:1810.04805_.
  * Gillick et al. (2018) Daniel Gillick, Alessandro Presta, and Gaurav Singh Tomar. 2018.  End-to-end retrieval in continuous space.  _arXiv preprint arXiv:1811.08008_.
  * Goodfellow et al. (2013) Ian Goodfellow, David Warde-Farley, Mehdi Mirza, Aaron Courville, and Yoshua Bengio. 2013.  Maxout networks.  In _International conference on machine learning_ , pages 1319–1327. PMLR.
  * Guu et al. (2015) Kelvin Guu, John Miller, and Percy Liang. 2015.  Traversing knowledge graphs in vector space.  _arXiv preprint arXiv:1506.01094_.
  * Hamilton et al. (2018) Will Hamilton, Payal Bajaj, Marinka Zitnik, Dan Jurafsky, and Jure Leskovec. 2018\.  Embedding logical queries on knowledge graphs.  _Advances in neural information processing systems_ , 31.
  * He et al. (2020) Kaiming He, Haoqi Fan, Yuxin Wu, Saining Xie, and Ross Girshick. 2020.  Momentum contrast for unsupervised visual representation learning.  In _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_ , pages 9729–9738.
  * Hu et al. (2022) Zhiwei Hu, Víctor Gutiérrez-Basulto, Zhiliang Xiang, Xiaoli Li, Ru Li, and Jeff Z Pan. 2022.  Type-aware embeddings for multi-hop reasoning over knowledge graphs.  _arXiv preprint arXiv:2205.00782_.
  * Kotnis et al. (2021) Bhushan Kotnis, Carolin Lawrence, and Mathias Niepert. 2021.  Answering complex queries in knowledge graphs with bidirectional sequence encoders.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 35, pages 4968–4977.
  * Lin et al. (2015) Yankai Lin, Zhiyuan Liu, Huanbo Luan, Maosong Sun, Siwei Rao, and Song Liu. 2015\.  Modeling relation paths for representation learning of knowledge bases.  _arXiv preprint arXiv:1506.00379_.
  * Liu et al. (2019) Yinhan Liu, Myle Ott, Naman Goyal, Jingfei Du, Mandar Joshi, Danqi Chen, Omer Levy, Mike Lewis, Luke Zettlemoyer, and Veselin Stoyanov. 2019.  Roberta: A robustly optimized bert pretraining approach.  _arXiv preprint arXiv:1907.11692_.
  * Ren et al. (2020) Hongyu Ren, Weihua Hu, and Jure Leskovec. 2020.  Query2box: Reasoning over knowledge graphs in vector space using box embeddings.  _arXiv preprint arXiv:2002.05969_.
  * Ren and Leskovec (2020) Hongyu Ren and Jure Leskovec. 2020.  Beta embeddings for multi-hop logical reasoning in knowledge graphs.  _Advances in Neural Information Processing Systems_ , 33:19716–19726.
  * Suchanek et al. (2007) Fabian M Suchanek, Gjergji Kasneci, and Gerhard Weikum. 2007.  Yago: a core of semantic knowledge.  In _Proceedings of the 16th international conference on World Wide Web_ , pages 697–706.
  * Teru et al. (2020) Komal Teru, Etienne Denis, and Will Hamilton. 2020.  Inductive relation prediction by subgraph reasoning.  In _International Conference on Machine Learning_ , pages 9448–9457. PMLR.
  * Toutanova and Chen (2015) Kristina Toutanova and Danqi Chen. 2015.  Observed versus latent features for knowledge base and text inference.  In _Proceedings of the 3rd workshop on continuous vector space models and their compositionality_ , pages 57–66.
  * Vrandečić and Krötzsch (2014) Denny Vrandečić and Markus Krötzsch. 2014.  Wikidata: a free collaborative knowledgebase.  _Communications of the ACM_ , 57(10):78–85.
  * Wang et al. (2021a) Bo Wang, Tao Shen, Guodong Long, Tianyi Zhou, Ying Wang, and Yi Chang. 2021a.  Structure-augmented text representation learning for efficient knowledge graph completion.  In _Proceedings of the Web Conference 2021_ , pages 1737–1748.
  * Wang et al. (2022) Liang Wang, Wei Zhao, Zhuoyu Wei, and Jingming Liu. 2022.  Simkgc: Simple contrastive knowledge graph completion with pre-trained language models.  _arXiv preprint arXiv:2203.02167_.
  * Wang et al. (2021b) Xiaozhi Wang, Tianyu Gao, Zhaocheng Zhu, Zhengyan Zhang, Zhiyuan Liu, Juanzi Li, and Jian Tang. 2021b.  Kepler: A unified model for knowledge embedding and pre-trained language representation.  _Transactions of the Association for Computational Linguistics_ , 9:176–194.
  * Wolf et al. (2019) Thomas Wolf, Lysandre Debut, Victor Sanh, Julien Chaumond, Clement Delangue, Anthony Moi, Pierric Cistac, Tim Rault, Rémi Louf, Morgan Funtowicz, et al. 2019.  Huggingface’s transformers: State-of-the-art natural language processing.  _arXiv preprint arXiv:1910.03771_.
  * Xiong et al. (2017) Wenhan Xiong, Thien Hoang, and William Yang Wang. 2017.  Deeppath: A reinforcement learning method for knowledge graph reasoning.  _arXiv preprint arXiv:1707.06690_.
  * Yao et al. (2019) Liang Yao, Chengsheng Mao, and Yuan Luo. 2019.  Kg-bert: Bert for knowledge graph completion.  _arXiv preprint arXiv:1909.03193_.
  * Zha et al. (2022) Hanwen Zha, Zhiyu Chen, and Xifeng Yan. 2022.  Inductive relation prediction by bert.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 36, pages 5923–5931.
  * Zhang et al. (2021) Zhanqiu Zhang, Jie Wang, Jiajun Chen, Shuiwang Ji, and Feng Wu. 2021.  Cone: Cone embeddings for multi-hop reasoning over knowledge graphs.  _Advances in Neural Information Processing Systems_ , 34:19172–19183.

##  Appendix A Structure Linearization Templates

In this part, we list the detailed linearization templates of different
structural query types in Table 7.

Query Type |  Linearization Template
---|---
[1p] |  [projection] [anchor] $t(e_{1})$ [relation] $t(r_{1})$
[2p] |  [projection] [anchor] $t(e_{1})$ [relation] $t(r_{1})$ [projection] [relation] $t(r_{2})$
[3p] |  [projection] [anchor] $t(e_{1})$ [relation] $t(r_{1})$ [projection] [relation] $t(r_{2})$ [projection] [relation] $t(r_{3})$
[2i] |  [intersection] [projection] [anchor] $t(e_{1})$ [relation] $t(r_{1})$ [projection] [anchor] $t(e_{2})$ [relation] $t(r_{2})$
[3i] |  [intersection] [projection] [anchor] $t(e_{1})$ [relation] $t(r_{1})$ [projection] [anchor] $t(e_{2})$ [relation] $t(r_{2})$ [projection] [anchor] $t(e_{3})$ [relation] $t(r_{3})$
[pi] |  [intersection] [projection] [anchor] $t(e_{1})$ [relation] $t(r_{1})$ [projection] [relation] $t(r_{2})$ [projection] [anchor] $t(e_{2})$ [relation] $t(r_{3})$
[ip] |  [intersection] [projection] [anchor] $t(e_{1})$ [relation] $t(r_{1})$ [projection] [relation] $t(r_{3})$ [projection] [anchor] $t(e_{2})$ [relation] $t(r_{2})$ [projection] [relation] $t(r_{3})$
[2u] |  [union] [projection] [anchor] $t(e_{1})$ [relation] $t(r_{1})$ [projection] [anchor] $t(e_{2})$ [relation] $t(r_{2})$
[up] |  [union] [projection] [anchor] $t(e_{1})$ [relation] $t(r_{1})$ [projection] [relation] $t(r_{3})$ [projection] [anchor] $t(e_{2})$ [relation] $t(r_{2})$ [projection] [relation] $t(r_{3})$

Table 7:  The linearization templates of each query structure.

Dataset | Training | Validation | Test
---|---|---|---
1p | others | 1p | others | 1p | others
FB15k | 273,710 | 273,710 | 59,097 | 8,000 | 67,016 | 8,000
FB15k-237 | 149,689 | 149,689 | 20,101 | 5,000 | 22,812 | 5,000
NELL995 | 107,982 | 107,982 | 16,927 | 4,000 | 17,034 | 4,000
FB15k-237-V2 | 9,964 | 9,964 | 1,738 | 2,000 | 791 | 1,000
NELL-V3 | 12,010 | 12,010 | 2,197 | 2,000 | 1,167 | 1,500

Table 8:  Statistics of logical queries in different types in each data split.
“others” means other eight query types including 2p, 3p, 2i, 3i, pi, ip, 2u
and up.

##  Appendix B Data Statistics

In Table 8, we summarize the statistics of logical queries in our experimented
datasets for both inductive and transductive settings.

##  Appendix C Constructed Query Types

We here introduce our generated more complicated query types involving more
relations and logical operations, which are used to illustrate the complicated
query structure generalizability.

$\displaystyle[4p]\ q=v_{?}$ $\displaystyle.\exists
v_{1},v_{2},v_{3}:r_{1}(e_{1},v_{1})\wedge r_{2}(v_{1},v_{2})$
$\displaystyle\wedge r_{3}(v_{2},v_{3})\wedge r_{4}(v_{3},v_{?})$
$\displaystyle[5p]\ q=v_{?}$ $\displaystyle.\exists
v_{1},v_{2},v_{3},v_{4}:r_{1}(e_{1},v_{1})\wedge r_{2}(v_{1},v_{2})$
$\displaystyle\wedge r_{3}(v_{2},v_{3})\wedge r_{4}(v_{3},v_{4})\wedge
r_{5}(v_{4},v_{?})$ $\displaystyle[3ip]\ q=v_{?}$ $\displaystyle.\exists
v:r_{1}(e_{1},v)\wedge r_{2}(e_{2},v)\wedge r_{3}(e_{3},v)$
$\displaystyle\wedge r_{4}(v,v_{?})$ $\displaystyle[i2p]\ q=v_{?}$
$\displaystyle.\exists v_{1},v_{2}:r_{1}(e_{1},v_{1})\wedge
r_{2}(e_{2},v_{1})$ $\displaystyle\wedge r_{3}(v_{1},v_{2})\wedge
r_{4}(v_{2},v_{?})$


# Lattice protein design using Bayesian learning

Tomoei Takahashi<EMAIL_ADDRESS>Graduate School of
Informatics, Nagoya University, Nagoya, 464-8601, Japan.  George Chikenji
Graduate School of Engineering, Nagoya University, Nagoya, 464-8603, Japan.
Kei Tokita  Graduate School of Informatics, Nagoya University, Nagoya,
464-8601, Japan.

###### Abstract

A novel protein design method using Bayesian learning is proposed in this
work. We consider a posterior probability of amino acid sequences by taking
into account water and assuming a prior of sequences. For some instances of a
target conformation of a two-dimensional (2D) lattice Hydrophobic-Polar (HP)
model, our method successfully finds an amino acid sequence for which the
target conformation has a unique ground state. However, the performance was
not as good for 3D lattice HP models compared with 2D models. Furthermore, we
find a strong linearity between the chemical potential of water and the number
of surface residues, thereby revealing the relationship between protein
structure and the effect of water molecules. The advantage of our method is
that it greatly reduces computation time, because it does not require long
calculations for the partition function corresponding to an exhaustive
conformational search. As our method uses a general form of Bayesian learning
and statistical mechanics and is not limited to lattice HP proteins, the
results presented here elucidate some heuristics used successfully in previous
protein design methods.

††preprint: APS/123-QED

##  I INTRODUCTION

Proteins have important roles in living systems. The complex three-dimensional
(3D) structure that determines the function of a protein represents an
equilibrium state determined by the amino acid sequence and the appropriate
physiological conditionsanfinsen1973principles . Protein
designcoluzza2017computational  requires determining the optimal amino acid
sequence that results in a given 3D structure as an equilibrium state. Protein
design, therefore, is the inverse problem of 3D structure prediction. In
addition, as amino acid sequences are determined by genomic information,
protein design can be used to explore the design principles of life. In recent
decades, many computational protein design methods have been proposed and
applied to drug designdahiyat1997novo ; kuhlman2003design ;
fleishman2011computational ; koga2012principles ; bale2016accurate ;
silva2019novo ; liu2006rosettadesign ; huang2011rosettaremodel . However,
there have been fewer theoretical studies of design principles based on
statistical mechanics, and few heuristics have been applied to the design of
real proteins.

In order to design an optimal sequence, the most reasonable statistical
mechanical procedure involves finding a sequence that minimizes the free
energy of a given target conformation. As part of this procedure, one needs to
carry out a folding simulation to check that the selected sequence folds into
the target conformation with high probability every time a candidate optimal
sequence is selected, for example, using the negative design methodjin2003novo
.

However, many other design methods stabilize only a target conformation,
although they could in principle lower the energy of other conformations.
Nevertheless, they successfully find an optimal sequence that minimizes the
energy of the target conformation. It is not clear why such design (inverse
problem) methods work without an exhaustive search (forward problem) for all
possible compact conformations. One of the main purposes of the present study
is to investigate the success of the inverse problem without solving the
forward problem, which is a highly non-trivial problem in statistical and
biological physics.

##  II MODEL AND METHOD

###  II.1 The lattice HP model

In order to address the problem described above, we take a statistical
mechanical approachcocco2018inverse  based on Bayesian learning, using a
coarse-grained protein model called the HP modellau1989lattice . We use a
lattice model in which every amino acid residue is located on a lattice site,
and a protein structure is represented by a self-avoiding walk on a 2D or 3D
lattice. Although a real protein has 20 types of amino acid, the HP model
includes only two types, hydrophobic (H) and polar (P). Here, we consider $N$
residues ${\bm{\sigma}}=\\{\sigma_{1},\sigma_{2},\ldots,\sigma_{N}|\forall
i,\sigma_{i}=\pm 1\\}$ on a lattice position
$\bm{r}=\\{r_{1},r_{2},\ldots,r_{N}\\}$, where $i=1,2,\ldots,N$ $\sigma_{i}=1$
indicates that the $i$-th residue is an H residue, and $\sigma_{i}=-1$
indicates that it is a P residue.

We assume the energy of the lattice protein is given by

$\displaystyle
E(\mbox{\boldmath$r$};{\bm{\sigma}})=\sum_{i<j}U(\sigma_{i},\sigma_{j})\Delta(r_{i}-r_{j}),$
(1)

where $U(\sigma_{i},\sigma_{j})$ denotes the interaction potential between the
monomers $i$ and $j$. We moreover assume the simplest functional form of $U$:
$U(1,1)=\epsilon_{1},U(1,-1)=U(-1,1)=\epsilon_{2},U(-1,-1)=\epsilon_{3}$,
using two types of interaction set,
$(\epsilon_{1},\epsilon_{2},\epsilon_{3})=(-1,0,0)$ and $(-2.3,-1,0)$. The
definition of the contact energy $\Delta(r_{i}-r_{j})$ is

$\displaystyle\Delta(r_{i}-r_{j})=\begin{cases}1&\text{ if $r_{i}$ and $r_{j}$
contact each other,}\\\ 0&\text{otherwise,}\end{cases}$ (2)

where contact between two residues is defined as the case where
$|r_{i}-r_{j}|=1$ but $|i-j|\neq 1$. In this model, therefore, the energy
given by Eq. (1) of denatured conformations is always higher than the energy
of compact conformations; thus, we only consider compact conformations in the
present study. Equilibrium statistical mechanics have been successfully
applied using the HP model, which is similar to the Ising model. For lattice
protein models of comparatively small size, several successful theoretical
studies have been reported.

###  II.2 Related works

The first and pioneering study of protein design using statistical mechanics
and the HP model was performed by Shakhnovich and Gutin
(SG)shakhnovich1993engineering . Subsequently, Kurosky and Deutsch proposed a
design criterion in which the solution of the design problem is a sequence
$\bm{\sigma}$ that maximizes the Boltzmann distribution
$p(\bm{R}|\bm{\sigma})$ as the conditional probabilitykurosky1995design ;
deutsch1996new , where $\bm{R}$ is a set of position vectors of the target
conformation. We hereafter call $p(\bm{R}|\bm{\sigma})$ the ‘target
probability,’ and the design criterion maximizing the target energy
probability is denoted the MTP criterion.

For the MTP criterion, a solution $\bm{\sigma}_{\rm MTP}$ is given by

$\displaystyle\bm{\sigma}_{\rm MTP}=\mathop{\rm
arg~{}max}\limits_{\bm{\sigma}}p(\bm{R}|\bm{\sigma})$ (3)

and

$\displaystyle p(\bm{R}|\bm{\sigma})=\frac{\exp\left[-\beta
E(\bm{R};{\bm{\sigma}})\right]}{Z_{\beta}({\bm{\sigma}})},$ (4) $\displaystyle
Z_{\beta}(\bm{\sigma})=\sum_{\bm{r}}\exp\left[-\beta
E(\bm{r};\bm{\sigma})\right],$ (5)

where $\beta$ is the inverse temperature of the system, $\bm{r}$ is any
compact conformation of the conformation space into which the sequence
$\bm{\sigma}$ can fold, and $Z_{\beta}({\bm{\sigma}})$ is a partition function
of the conformation space. In terms of statistical learning, the MTP criterion
is the maximum likelihood estimation for the canonical distribution [Eq. (4)].
The MTP criterion, however, includes $Z_{\beta}(\bm{\sigma})$. Hence, to
obtain $\bm{\sigma}_{\rm MTP}$, one must carry out an exhaustive
conformational search every time a candidate $\bm{\sigma}_{\rm MTP}$ is found.
For such conformational searches, very fast and accurate methods are required.
Currently available methods include generalized ensemble Monte Carlo
methodsiba1998simulation ; chikenji1999multi ; ShiraiKikuchi2013 . Even these
methods, however, cannot provide reasonable results for longer chains (more
than 100 residues). Thus, design of large and realistic protein models remains
impossible, and successful statistical mechanical protein design using the MTP
criterion has been reported only for comparatively small lattice proteins and
modelskurosky1995design ; deutsch1996new ; seno1996optimal ;
micheletti1998protein ; irback1998monte ; irback1999design ; iba1998design ;
tokita2000dynamical ; rossi2000anovel ; betancourt2002protein ; zou2003using ;
wang2004 ; jiao2006protein ; kleinman2006maximum ; coluzza2003designing ;
coluzza2011coarse .

To overcome the above difficulty, Kurosky and Deutschdeutsch1996new  carried
out high-temperature expansion of the free energy
$F_{\beta}(\bm{\sigma})=-(1/\beta)\ln{Z_{\beta}(\bm{\sigma})}$ and minimized
$F_{\beta}(\bm{\sigma})$ using simulated annealing for a 2D HP model with
$N=10$–$18$. A design method using simulated annealing for both sequence and
conformation space was proposed by Seno et al.seno1996optimal . using a 2D
lattice HP model with $N=12$ and $16$. The multi-sequence Monte Carlo method
proposed by Irbäck et al.irback1998monte ; irback1999design  is an efficient
procedure that obtains an optimal sequence excluding bad sequences with low
target probability using fluctuation of sequences; this method was used to
design a comparatively large 2D HP model ($N=32$ and $50$) and 3D off-lattice
HP model ($N=20$). The ‘design equation’ method of Iba et al. was the first
application of Boltzmann machine learning to protein design and was used to
obtain a correct sequence for several 3D $3\times 3\times 3$ cubic
conformationsiba1998design ; tokita2000dynamical . There are some other design
methods with a similar policy to the design equation method: minimization of
the relative entropywang2004 ; jiao2006protein  and gradient
descentkleinman2006maximum  to optimize the sequence. Coluzza et al. proposed
a method using both energy minimization and free energy minimization for some
targetscoluzza2003designing ; coluzza2011coarse .

In recent years, many successful applications of deep learning have been
reported in various engineering and scientific fields, including protein
folding and drug designAlphaFold ; li2014direct ; wang2018computational ;
greener2018design ; merk2018novo ; gupta2018generative ; o2018spin2 ;
button2019automated . Deep learning for protein design involves learning the
relations between a protein conformation and an amino acid sequence using big
data from the Protein Data Bankbernstein1977protein . Although deep neural
networks have been used to successfully predict an optimal sequence for a
target conformationwang2018computational ; button2019automated , we still do
not theoretically understand the design principles of machine learning in this
context. In the present study, therefore, we apply Bayesian learning to
protein design to explore these design principles.

###  II.3 Design method by Bayesian learning

In Bayesian learning, one assumes that data $D$ are generated by conditional
probability $p(D|\theta)$ under the parameter $\theta$. By Bayes’ theorem, the
posterior distribution $p(\theta|D)$ is given by

$\displaystyle p(\theta|D)=\frac{p(D|\theta)p(\theta)}{\int
p(D|\theta)p(\theta)d\theta},$ (6)

where $p(\theta)$ is the prior distribution of $\theta$ if $\theta$ is a
continuous random variable. One estimates $\theta$ and predicts unobserved
data from the posterior given by Eq. (6). The basic procedure of Bayesian
learning is as follows: starting from a highly arbitrary prior $p(\theta)$,
the posterior $p(\theta|D)$ is repeatedly updated with the highly objective
likelihood function $p(D|\theta)$ using Eq. (6); thus we can finally obtain a
precise value of $\theta$.

Here we apply the above procedure to protein design. Let the appearance
probability of the given target conformation $\bm{R}$ for data $D$ be the
target probability $p(\bm{R}|\bm{\sigma})$, and let the prior of the HP
sequence $\bm{\sigma}$ be the parameter $\theta$, $p(\bm{\sigma})$. Then, the
posterior of the sequence $p(\bm{\sigma}|\bm{R})$ is given by

$\displaystyle
p({\bm{\sigma}}|\bm{R})=\frac{p(\bm{R}|\bm{\sigma})p(\bm{\sigma})}{\sum_{\bm{\sigma}}p(\bm{R}|\bm{\sigma})p(\bm{\sigma})},$
(7)

where $\sum_{\bm{\sigma}}$ denotes summation over all sequences.

Furthermore, we focus on the fact that a protein acquires its native state as
an equilibrium with water molecules in a living cell.

We therefore consider the following target probability as a grand canonical
distribution:

$\displaystyle p(\bm{R},N_{\rm
w}|\bm{\sigma})=\frac{\exp\left[-\beta\left(E(\bm{R};\bm{\sigma})-\mu N_{\rm
w}\right)\right]}{\Xi_{\beta,\mu}(\bm{\sigma})},$ (8)
$\displaystyle\Xi_{\beta,\mu}(\bm{\sigma})=\sum_{N_{\rm
w}=0}^{\infty}\sum_{\bm{r}}\exp\left[-\beta\left(E(\bm{r};\bm{\sigma})-\mu
N_{\rm w}\right)\right].$ (9)

The definition of the energy $E(\bm{R};\bm{\sigma})$ of the target $\bm{R}$
for a given sequence $\bm{\sigma}$ is given by Eq. (1), and $\mu$ is the
chemical potential of water. We assume that one water molecule combines with
one P residue, hence $N_{\rm w}=N_{\rm P}(\bm{\sigma})$, where $N_{\rm
P}(\bm{\sigma})$ is the number of all P residues. Therefore, Eqs. (8) and (9)
denote the canonical distribution of the Hamiltonian
$E(\bm{R};\bm{\sigma})-\mu N_{\rm P}(\bm{\sigma})$ with external field $\mu$.
Consequently, we rewrite Eqs. (8) and (9) as follows:

$\displaystyle
p(\bm{R}|\bm{\sigma})=\frac{\exp\left[-\beta\left(E(\bm{R};\bm{\sigma})-\mu
N_{\rm P}(\bm{\sigma})\right)\right]}{\Xi_{\beta,\mu}(\bm{\sigma})},$ (10)
$\displaystyle\Xi_{\beta,\mu}(\bm{\sigma})=\sum_{\bm{r}}\exp\left[-\beta\left(E(\bm{r};\bm{\sigma})-\mu
N_{\rm P}(\bm{\sigma})\right)\right],$ (11)

where $N_{\rm P}(\bm{\sigma})$ is obtained by the conditions $N=N_{\rm
H}+N_{\rm P}$ and $\sum_{i}\sigma_{i}=N_{\rm H}-N_{\rm P}$ as

$\displaystyle N_{\rm P}=\frac{1}{2}(N-\sum_{i}\sigma_{i}).$ (12)

In order to obtain an optimal sequence $\bm{\sigma}$ that maximizes
$p(\bm{R}|\bm{\sigma})$, we have to repeat the exhaustive conformational
search in Eq. (11) for each trial sequence $\bm{\sigma}$. In general, this
calculation takes an enormous amount of time.

One of the new ideas of the present study is that the prior distribution
$p(\bm{\sigma})$ is given by

$\displaystyle
p(\bm{\sigma})=\frac{\Xi_{\beta,\mu}(\bm{\sigma})}{\Xi_{\beta,\mu}},$ (13)
$\displaystyle\Xi_{\beta,\mu}=\sum_{\bm{\sigma}}\sum_{\bm{r}}\exp\left[-\beta\left(E(\bm{r};\bm{\sigma})-\mu
N_{{\rm P}}(\bm{\sigma})\right)\right],$ (14)

where $\Xi_{\beta,\mu}$ is the partition function of both conformation and
sequence space. The expression of the prior Eq. (13) is based on the following
hypothesis: as a result of evolution, the probability $p(\bm{\sigma})$ is
proportional to its partition function $\Xi_{\beta,\mu}(\bm{\sigma})$ so that
the free energy
$F_{\beta,\mu}(\bm{\sigma}):=-(1/\beta)\ln\Xi_{\beta,\mu}(\bm{\sigma})$ takes
a minimum. Note that Kurosky and Deutschkurosky1995design  assumed the equal
_a priori_ weights $p(\bm{\sigma})=1/N_{\rm s}$, where $N_{\rm s}$ is the
number of all HP sequences. By contrast, our method considers the above
postulation for the weight of the appearance of sequences, which is reasonable
from the viewpoint of thermodynamics and protein evolution. We obtain
posterior distribution $p(\bm{\sigma}|\bm{R})$ by substituting Eqs. (13) and
(10) into Eq. (7) and by canceling $\Xi_{\beta,\mu}(\bm{\sigma})$ out as
follows:

$\displaystyle
p(\bm{\sigma}|\bm{R})=\frac{\exp\left[-\beta\left(E(\bm{R};\bm{\sigma})-\mu
N_{{\rm P}}(\bm{\sigma})\right)\right]}{\Xi_{\beta,\mu}(\bm{R})},$ (15)
$\displaystyle\Xi_{\beta,\mu}(\bm{R})=\sum_{\bm{\sigma}}\exp\left[-\beta\left(E(\bm{R};\bm{\sigma})-\mu
N_{{\rm P}}(\bm{\sigma})\right)\right],$ (16)

where $\Xi_{\beta,\mu}(\bm{R})$ is a partition function of the sequence space
corresponding to the given target conformation $\bm{R}$. An important point
regarding Eq. (15) is that it no longer includes
$\Xi_{\beta,\mu}(\bm{\sigma})$. Thus, we can efficiently obtain an optimal
sequence using Eq. (15) simply by summation over $\bm{\sigma}$. This design
method without $\Xi_{\beta,\mu}(\bm{\sigma})$ is essentially same as a
procedure used previously to obtain more realistic protein
designsdahiyat1997novo ; kuhlman2003design ; fleishman2011computational ;
koga2012principles ; bale2016accurate ; silva2019novo ; liu2006rosettadesign ;
huang2011rosettaremodel .

The final step is to obtain an optimal HP sequence using Eq. (15).
Nevertheless, the exact calculation of Eq. (16) is difficult if the number of
residues $N$ is large. We thus utilize one of the simplest Markov-chain Monte
Carlo (MCMC) methods, Gibbs sampling. In this method, the sampling probability
of $\sigma_{i}$ of each Monte Carlo step (MCS) is a conditional probability of
$\sigma_{i}$ given other random variables. We thus obtain the following
sampling probability by substituting Eqs. (1) and (12) into Eq. (15).
Accordingly, the sampling probability of an H residue ($\sigma_{i}=1$) or P
residue ($\sigma_{i}=-1$) is given by

$\displaystyle p(\sigma_{i}=\pm 1|\bm{R};\bm{\sigma}_{\backslash
i})=\frac{1}{1+\mathrm{e}^{\pm{\beta}\\{\Delta
E_{i}(\bm{R};\bm{\sigma})+\mu\\}}},$ (17)

where $\bm{\sigma}_{\backslash
i}:=\\{\sigma_{1},\dots,\sigma_{i-1},\sigma_{i+1},\ldots,\sigma_{N}\\}$, a
vector of all random variables of residues except for the $i$-th residue
$\sigma_{i}$, and the double signs correspond. Let $\Delta
E_{i}(\bm{R};\bm{\sigma}):=\sum_{j\in
n(i)}\left(U(1,\sigma_{j})-U(-1,\sigma_{j})\right)$, where $n(i)$ denotes the
set of sites $j$ that are the nearest neighbors of $i$-th site except for
those along the chain ($j\neq i-1,i+1$). The random variables
$\bm{\sigma}_{\backslash i}$ have fixed realizations in the denominator and
the numerator of the right-hand side of Eq. (15) at every MCS. Thus, the
random variables that interact with the $i$-th residue $\sigma_{i}$ remain
only on the right-hand side of Eq. (17), because those fixed realizations,
except for the residues that interact with $\sigma_{i}$, are canceled out in
Eq. (15). We decide whether each residue is H or P using the expectation
$\left<\sigma_{i}\right>$; that is, $\sigma_{i}$ is H if
$\left<\sigma_{i}\right>>0$ and P otherwise. We also take the number of MCSs
until the estimated value does not change and let the burn-in be the leading
$1/5$ of all MCSs. In this study, the inverse temperature is set to $\beta=10$
for all conformations of all lattice models. On the other hand, we
heuristically set the chemical potential $\mu$ in order to design a unique
ground state by repeating the design experiment many times. The necessary and
sufficient condition for successful design is that the energy given by Eq. (1)
of the target conformation and the sequence designed corresponds to a unique
ground state of all possible compact conformations.

The formulation of our design method as described so far can be derived by
assuming a joint distribution given by

$\displaystyle
p(\bm{r};\bm{\sigma})=\frac{\exp\left[-\beta\left(E(\bm{r};\bm{\sigma})-\mu
N_{{\rm P}}(\bm{\sigma})\right)\right]}{\Xi_{\beta,\mu}},$ (18)

where $\Xi_{\beta,\mu}$ is given by Eq. (14). One can derive the prior [Eq.
(13)] and the likelihood function (the target probability) given by Eq. (10)
by a marginalization $p(\bm{\sigma})=\sum_{\bm{r}}p(\bm{r};\bm{\sigma})$ and a
relation
$p(\bm{R}|\bm{\sigma})=\left(p(\bm{r};\bm{\sigma})/p(\bm{\sigma})\right)|_{\bm{r}=\bm{R}}$,
respectively. Thus, the hypothesis [Eq. (13)] is included in the joint
distribution [Eq. (18)].

##  III RESULT

###  III.1 Enumerable conformations

First, we tested our design method with comparatively small lattice protein
models, for which all compact conformations were enumerable. We designed 2D
$N=3\times 3$, $3\times 4$, $4\times 4$, $5\times 5$, and $6\times 6$ lattice
models, and 3D $N=2\times 2\times 3$ and $3\times 3\times 3$ lattice models.
The numbers of all conformations $N_{\rm c}$ and the numbers of all HP
sequences $N_{\rm s}$ of these lattice models are shown in Table 1.

The number of conformations $N_{\rm c}$ is the number of all compact self-
avoiding walks from which all kind of rotational, reflection, and head-tail
symmetrical conformations have been eliminated. The total number of
conformations of these lattice models is enumerable; hence, one can confirm
whether or not the designed sequence folds into the target conformation as a
unique ground state. Note that not every conformation always has a solution to
the design problem. The number of sequences that fold into the target
conformation as a unique ground state differs among target conformations and

Table 1: Number of conformations and HP sequences Size | $N_{\rm c}$ | $N_{\rm s}$ | Conformations designed
---|---|---|---
$3\times 3$ | 8 | 512 | All
$3\times 4$ | 27 | 4096 | All
$4\times 4$ | 62 | 65536 | All
$5\times 5$ | 1075 | 33554432 | Random 100
$6\times 6$ | 52667 | 68719476736 | Random 100 and MHDC
$2\times 2\times 3$ | 69 | 4096 | All
$3\times 3\times 3$ | 103346 | 134217728 | Random 100 and MHDC

is called designability. In general, a conformation with higher designability
is easier to design. This is because high designability means a large solution
space in sequence space.

Designability is a significant quantity that relates to the thermodynamic
stability of proteins; however, we do not address issues of designability in
depth here. In order to calculate the exact success rate ($SR$) of the overall
conformation, one needs to select designable target conformations with
designability greater than zero; however, to enumerate the designabilities of
each conformation, one would need to enumerate the energy of every combination
of conformations and sequences. This would require vast computation time for
models with comparatively large size, such as the $5\times 5$, $6\times 6$,
and $3\times 3\times 3$ models (Table 1), even though they are compact.
Therefore, in this study, we carried out the enumeration of designabilities
only for the $3\times 3$, $3\times 4$, $4\times 4$, and $2\times 2\times 3$
lattice models.

For the models with $N=5\times 5$, $6\times 6$, and $3\times 3\times 3$, the
number of conformations was too large. Thus, we randomly chose 100 target
conformations and determined the $SR$, that is, the number of successfully
designed conformations (Table 2). For the models with $N=6\times 6$ and
$3\times 3\times 3$, we moreover identified the most highly designable
conformation (MHDC) (Figs. 1, 2, and 3), in which designabilities were exactly
enumeratedli1996emergence , to test whether our method could be used to design
the easiest instance.

The results of the application of our method are summarized in Table 2. All
designed sequences were classified into three types: good, medium, and bad
sequences. The good sequences had the target conformation as a unique ground
state, medium sequences had the target conformation as one of the degenerated
ground states, and bad sequences had ground state conformation(s) that did not
include the target conformation. In the table, $SR$, $N_{\rm c}^{(\rm g)}$,
$N_{\rm c}^{(\rm m)}$, and $N_{\rm c}^{(\rm b)}$, denote the percentage of
good sequences and the number of conformations that were designed with good,
medium, and bad sequences, respectively. We also calculated the average
degeneracy, $d_{\rm av}$, for all $N_{\rm c}^{(\rm g)}+N_{\rm c}^{(\rm m)}$
ground states. We repeated the calculations with various values of $\mu$ and
obtained the optimal value $\mu^{*}$ that gave the maximum success rate. The
values of $\mu^{*}$ for each lattice size are listed in Table 2. The values of
the energy parameters are also listed in Table 2. The energy parameters
$(\epsilon_{1},\epsilon_{2},\epsilon_{3})=(-2.3,-1,\,0)$ were also used for a
$3\times 3\times 3$ lattice model in previous workli1996emergence  in order to
avoid the degeneracy of ground states. We used the same energy parameters for
3D lattice models for the same reason. The total MCSs were set to $10^{5}$ for
all target conformations listed in Table 2. The small $3\times 3$ and $2\times
2\times 3$ lattices included several non-designable conformations; we excluded
such conformations when calculating $SR$.

According to the results shown in Table 2, success rates were relatively high
for small 2D HP models, but they decreased as $N$ increased. The average
degeneracy $d_{{\rm av}}$ was low for 2D models. By contrast, the success rate
for 3D models

Table 2: Design results and the optimal chemical potential $\mu^{*}$ Size | $N_{\rm c}^{(\rm g)}$ | $N_{\rm c}^{(\rm m)}$ | $N_{\rm c}^{(\rm b)}$ | $SR$ | $d_{{\rm av}}$ | $\mu^{*}$ | $(\epsilon_{1},\epsilon_{2},\epsilon_{3})$ |
---|---|---|---|---|---|---|---|---
$3\times 3$ | $4$ | $0$ | $0$ | 100 | 1 | $0.55$ |  |
$3\times 4$ | $20$ | 7 | 0 | $74.1$ | 1.33 | 0.6 |  |
$4\times 4$ | 46 | 16 | 0 | 74.2 | 1.32 | 0.62 | $(-1,0,0)$ |
$5\times 5$ | 68 | 32 | 0 | 68 | 1.48 | $0.74$ |  |
$6\times 6$ | 63 | 37 | 0 | $63$ | 1.58 | 0.8 |  |
$2\times 2\times 3$ | $17$ | 30 | 1 | 35.4 | 2.94 | 1.7 | $(-2.3,-1,0)$ |
$3\times 3\times 3$ | 8 | 80 | 12 | 8 | 10.67 | 2.33 |

was low compared with that of 2D models. For $2\times 2\times 3$, $d_{{\rm
av}}$ was low, but for $3\times 3\times 3$, it was comparatively high. Thus,
designed sequences did not appear to be likely to fold into the target
conformations for the $3\times 3\times 3$ cubic lattice. In addition,
$\mu^{*}$ increased as the number of residues increased for both the 2D and 3D
lattices.

Figure 1: Designed sequence of the MHDC of $6\times 6$ HP model with
$(\epsilon_{1},\epsilon_{2},\epsilon_{3})=(-1,0,0)$, $\beta=10$, and
$\mu^{*}=0.8$. The white and black balls denote H and P residues, respectively
(the same applies in the following figures. Figure 2:  Designed sequence of
the MHDC $3\times 3\times 3$ HP model with
$(\epsilon_{1},\epsilon_{2},\epsilon_{3})=(-2.3,-1,0)$, $\beta=10$, and
$\mu^{*}=2.33$. Figure 3:  Designed sequence of the MHDC of $3\times 3\times
3$ HP model with $(\epsilon_{1},\epsilon_{2},\epsilon_{3})=(-1,0,0)$,
$\beta=10$, and $\mu^{*}=1.0$.

Note that we did not enumerate designabilities of all conformations for the
$5\times 5$, $6\times 6$, and $3\times 3\times 3$ models; hence, there may
have been non-designable conformations among the 100 randomly chosen
conformations. However, it is likely that this was not the case for the
$5\times 5$ and $6\times 6$ models, because the smaller $3\times 4$ HP model
did not lead to any non-designable conformation. On the other hand, the
fraction of non-designable conformations out of all conformations for the
$2\times 2\times 3$ model was 21/69; the fraction for the $3\times 3\times 3$
model is expected to be less than that because the fraction decreased as the
size increased in the 2D cases. Thus, there may have been a considerable
number of non-designable conformations among the randomly chosen 100
conformations for the $3\times 3\times 3$ model; hence, the real success rate
of the $3\times 3\times 3$ model increased when non-designable conformations
were excluded.

Concerning the MHDC of the $6\times 6$ and $3\times 3\times 3$ HP models, we
obtained a good sequence (Figs. 1, 2, and 3). This is the first example of
design of a $6\times 6$ MHDC without enumerating all HP
sequencesli1996emergence . For the MHDC of the $3\times 3\times 3$ HP model,
we successfully designed a good sequence for the energy parameters
$(\epsilon_{1},\epsilon_{2},\epsilon_{3})=(-2.3,-1,\,0)$ (Fig. 2) and (-1, 0,
0) (Fig. 3). We executed $10^{4}$ MCSs for these three cases. The results
obtained here demonstrate the features of general globular proteins, with H
residues on the inside of the protein and P residues on the surface exposed to
the surrounding water molecules. We observed four residues (surrounded by
dotted black circles in Figs. 2 and 3) that were different from each other,
possibly owing to the presence or absence of H-P (P-H) contact energies.

The larger $\mu^{*}$ of the MHDC of $3\times 3\times 3$ HP model with
$(\epsilon_{1},\epsilon_{2},\epsilon_{3})=(-2.3,-1,\,0)$ compared with the
case of $(\epsilon_{1},\epsilon_{2},\epsilon_{3})=(-1,0,0)$ could have been
due to the lower H-H interaction $\epsilon_{1}=-2.3$, leading to a greater
increase in the number of H-residues than in the case of
$(\epsilon_{1},\epsilon_{2},\epsilon_{3})=(-1,0,0)$. Therefore, one needs to
increase $\mu^{*}$ in order for the surface residues to be P residues.

###  III.2 Large 2D conformations

Here, we chose 2D HP models with comparatively large size ($N=32,50$) models
studied by Irbäck et al.irback1998monte ; irback1999design . This confirmed
that the designed sequence was likely to fold into the target conformation
with simulated tempering. For the model with $N=32$ (respectively 50), the
parameters were set to $\mu^{*}=0.7$ (0.85) and the MCSs were $10^{4}$
($10^{5}$). The energy parameters were set to
$(\epsilon_{1},\epsilon_{2},\epsilon_{3})=(-1,0,0)$ in both cases. The
simulation was executed by a normal PC with 1.2 GHz dual-core Intel Core m3
and 8 GB memory, and the calculation time was approximately 0.5–1 s (11–12 s)
for $N=32$ (50). Thus, our method ran faster than those used in previous
studies. As a result, we successfully designed the same sequences reported by
Irbäck et al. (Figs. 4 and 5). Our method also demonstrates the features of
globular proteins.

Figure 4:  Designed sequence of the $N=32$ 2D HP model with
$(\epsilon_{1},\epsilon_{2},\epsilon_{3})=(-1,0,0)$, $\beta=10$, and
$\mu^{*}=0.7$.  Figure 5:  Designed sequence of the $N=50$ 2D HP model with
$(\epsilon_{1},\epsilon_{2},\epsilon_{3})=(-1,0,0)$, $\beta=10$, and
$\mu^{*}=0.85$.  Figure 6: Relation between $\mu^{*}$ and number of surface
residues, $N_{\rm sur}$.

###  III.3 Optimal $\mu^{*}$ and number of surface residues

We represent the relation between the optimal $\mu^{*}$ and the number of
surface residues $N_{\rm sur}$ in Fig. 6. We show only the results for
$(\epsilon_{1},\epsilon_{2},\epsilon_{3})=(-1,0,0)$ because the optimal
$\mu^{*}$ varies depending on the energy parameters for a given conformation.
We therefore plotted the results for all 2D models and the $3\times 3\times 3$
MHDC model with $(\epsilon_{1},\epsilon_{2},\epsilon_{3})=(-1,0,0)$ The
residues that were bent 90 degrees inward (indicated by a dashed black circle
in Figs. 4 and 5) were not counted for $N_{\rm sur}$ because a water molecule
is unlikely to combine with such residues (see Fig. 6).

We observed noticeable linearity between $\mu^{*}$ and $N_{\rm sur}$. The
outlier $(\mu^{*},N_{\rm sur})=(0.70,20)$ was obtained in the 2D $N=32$ case
(Fig. 4), in which the target conformation was not fully compact and the
number of surface residues was much larger than those of other target
conformations tested. According to these results, the optimal $\mu^{*}$ can be
estimated by the number of surface residues of a target conformation.

###  III.4 Probability of a P residue

Finally, in order to clarify why 3D models performed less well than 2D models,
we consider the probability $P_{\rm P}$ that a residue is P for all residues
of the $3\times 3\times 3$ and $6\times 6$ MHDC models (Figs. 7 and 8). We use
$p(\sigma_{i}=-1|\bm{R};\bm{\sigma}_{\backslash i})$ in Eq. (17) as $P_{\rm
P}$.

Figure 7: $P_{\rm P}$ of each residue of $3\times 3\times 3$ MHDC with
$(\epsilon_{1},\epsilon_{2},\epsilon_{3})=(-2.3,-1,\,0)$ and $\mu=2.33$ (Fig.
2). The residue number starts from the center residue of the front side of
Fig. 2. Figure 8: $P_{\rm P}$ of each residue of $6\times 6$ MHDC with
$(\epsilon_{1},\epsilon_{2},\epsilon_{3})=(-1,0,0)$ and $\mu=0.8$ (Fig. 1).
The residue number starts from the bottom left residue of Fig. 1.

Each $P_{\rm P}$ in Figs. 7 and 8 is the average of
$1/\left(1+\exp\left[-{\beta}\\{\Delta
E_{i}(\bm{R};\bm{\sigma})+\mu\\}\right]\right)$ over the last $1/5$ MCSs in
both cases. In the case of the 3D $3\times 3\times 3$ lattice, the values of
$P_{\rm P}$ for residues 2, 16, 22, and 24 were not very high. These residues
were located in the center of each cube side. By contrast, all $P_{\rm P}$
values greater than 0.5 were almost equal to 1 in the case of the 2D $6\times
6$ lattice. Accordingly, it can be seen that the clear division of all $P_{\rm
P}$ values into $1$ or $0$ is an index of successful design. Thus, the 3D
models are difficult instances for our design method.

##  IV DISCUSSION

Our method is similar to the SG method proposed by Shakhnovich and
Gutinshakhnovich1993engineering , which did not include $Z(\bm{\sigma})$. The
difference between the SG method and ours is the minimization function. The SG
method minimizes $E(\bm{R};\bm{\sigma})$ directly, keeping $N_{\rm
P}(\bm{\sigma})$ a constant value determined _a priori_ , but our method
minimizes $E(\bm{R};\bm{\sigma})-\mu N_{\rm P}(\bm{\sigma})$. Therefore, our
method can minimize $E(\bm{R};\bm{\sigma})$, maintaining the general features
of globular proteins, that is, H residues on the inside and P residues on the
surface exposed to the surrounding water molecules. Thus, one can minimize
$E(\bm{R};\bm{\sigma})$ while reducing the diversity of conformations into
which a designed HP sequence can fold by minimizing $E(\bm{R};\bm{\sigma})-\mu
N_{\rm P}(\bm{\sigma})$. This corresponds, in a sense, to negative
designjin2003novo .

On the other hand, as discussed above, our method failed in the cases of 3D HP
models and comparatively large compact 2D HP models because it failed to
reduce the diversity of the foldable conformations of the designed sequence in
these cases. The diversity—or, more simply, the total number of self-avoiding
walks—increases for 3D lattices compared with 2D lattices. Such diversity is
expected to increase as $N$ increases even in the case of 2D lattices. Thus,
the success rates decrease in these cases.

In addition, the numbers of core residues in the 3D models used in this study
were low, e.g., the 3D $2\times 2\times 3$ model had no core residue and the
3D $3\times 3\times 3$ model had only one; hence, it was difficult to design
globular protein-like conformations using these models. Given that our design
method finds an optimal sequence by controlling $\mu$, such small numbers of
cores may explain the lower performance of our method in the case of the 3D
models. If our design method works in a given instance, the posterior [Eq.
(15)] should show a sharp peak at the optimal $\mu^{*}$. This is equivalent to
the case where the $P_{\rm P}$ of every residue is almost equal to 1 or 0. For
the 3D $3\times 3\times 3$ lattice, however, there were several comparatively
low $P_{\rm P}$ values (close to 0.5), even in the case of the highly
designable conformation. Thus, our design method was not appropriate for those
conformations.

The greatest advantage of our method is that it skips the exhaustive
calculation of $Z(\bm{\sigma})$ by assuming the prior distribution given in
Eq. (13). As already stated, the form of the prior means that the lower the
free energy
$F_{\beta,\mu}(\bm{\sigma}):=-(1/\beta)\ln\Xi_{\beta,\mu}(\bm{\sigma})$, the
higher the prior distribution. The prior [Eq. (13)] states the hypothesis that
sequences rich in P residues are, in general, more likely to be evolutionarily
selected than sequences with unique stable conformations. This hypothesis is
consistent with recent findings that organisms have many intrinsically
disordered proteinswright1999intrinsically ; uversky2000natively ;
dunker2001intrinsically ; tompa2002intrinsically ; these proteins do not have
unique native conformations and are composed of sequences rich in P residues.
Recent workuversky2017intrinsically  has shown that such proteins form
‘droplets (that function as membraneless organelles)’ and have various
biologically important roles (e.g., spatiotemporal regulation of gene
expression, signaling, and stress response). This indicates that organisms
make good use of the physical property given by Eq. (13).

In addition to the biological validity of the prior, the fact that it enabled
fast protein design without the calculation of $Z(\bm{\sigma})$ is significant
because it suggests that Eq. (13) is not a unique solution for protein design
without exhaustive calculation of $Z(\bm{\sigma})$. As all information about
the thermodynamic profile of a protein is evolutionarily embedded solely in
the sequence, it is in principle possible to search for a sequence that
stabilizes a given target conformation if the code connecting $\bm{\sigma}$
and the thermodynamic profile is broken. The prior given by Eq. (13) may be
one such code.

##  V CONCLUSION

The simple conclusion from these results is that it is possible to design many
conformations without an exhaustive conformational search by taking the water
effect into account. This approach is more successful with 2D HP models than
with 3D models; however, our method is expected to correctly design 3D target
conformations given a sufficiently high designability of the target
conformation. This is consistent with conventional protein design software,
e.g., Rosetta. Hence, our method based on statistical mechanics may enable
future studies on more realistic protein design.

Future work could consider an additional parameter reflecting the specific
topological information of a target conformation. In addition, setting
different numbers of water molecules to combine with each P residue would help
to more closely model realistic globular proteins. As Bayesian learning is
simple and flexible, such modifications could be easily implemented.

## ACKNOWLEDGMENTS

The authors are grateful to M. Ota, Nagoya University, for illuminating
discussions. This work was supported by KAKENHI Nos. 19H03166 (G.C.) and
19K03650 (K.T.).

## References

  * (1) C. B. Anfinsen, Science 181, 223 (1973).
  * (2) I. Coluzza, J. Phys. Condens. Matt. 29, 143001 (2017).
  * (3) B. I. Dahiyat and S. L. Mayo, Science 278, 82 (1997).
  * (4) B. Kuhlman, G. Dantas, G. C. Ireton, G. Varani, B. L. Stoddard, and D. Baker, Science 302, 1364 (2003).
  * (5) S. J. Fleishman, T. A. Whitehead, D. C. Ekiert, C. Dreyfus, J. E. Corn, E.-M. Strauch, I. A. Wilson, and D. Baker, Science 332, 816 (2011).
  * (6) N. Koga, R. Tatsumi-Koga, G. Liu, R. Xiao, T. B. Acton, G. T. Montelione, and D. Baker, Nature 491, 222 (2012).
  * (7) J. B. Bale, S. Gonen, Y. Liu, W. Sheffler, D. Ellis, C. Thomas, D. Cascio, T. O. Yeates, T. Gonen, N. P. King, and D. Baker, Science 353, 389 (2016).
  * (8) D.-A. Silva, S. Yu, U. Y. Ulge, J. B. Spangler, K. M. Jude, C. Labão-Almeida, L. R. Ali, A. Quijano-Rubio, M. Ruterbusch, I. Leung, T. Biary, S. J. Crowley, E. Marcos, C. D. Walkey, B. D. Weitzner, F. Pardo-Avila, J. Castellanos, L. Carter, L. Stewart, S. R. Riddell, M. Pepper, G. J. L. Bernardes, M. Dougan, K. C. Garcia, and D. Baker, Nature 565, 186 (2019).
  * (9) Y. Liu and B. Kuhlman, Nucleic Acids Res. 34, W235 (2006).
  * (10) P.-S. Huang, Y.-E. A. Ban, F. Richter, I. Andre, R. Vernon, W. R. Schief, and D. Baker, PLoS ONE 6, e24109 (2011).
  * (11) W. Jin, O. Kambara, H. Sasakawa, A. Tamura, and S. Takada, Structure 11, 581 (2003).
  * (12) S. Cocco, C. Feinauer, M. Figliuzzi, R. Monasson, and M. Weigt. Rep. Prog. Phys. 81, 032601 (2018).
  * (13) K. F. Lau and K. A. Dill, Macromolecules 22, (1989).
  * (14) E. I. Shakhnovich and A. M. Gutin, Proc. Natl. Acad. Sci. U.S.A. 90, 719 (1993).
  * (15) T. Kurosky and J. M. Deutsch, J. Phys. A 28, L387 (1995).
  * (16) J. M. Deutsch and T. Kurosky, Phys. Rev. Lett. 76, 323 (1996).
  * (17) Y. Iba, G. Chikenji, and M. Kikuchi, J. Phys. Soc. Jpn. 67, 3327 (1998).
  * (18) G. Chikenji, M. Kikuchi, and Y. Iba, Phys. Rev. Lett. 83, 1886 (1999).
  * (19) N. C. Shirai and M. Kikuchi, J. Chem. Phys. 139, 225103-1 (2013).
  * (20) F. Seno, M. Vendruscolo, A. Maritan, and J. R. Banavar, Phys. Rev. Lett. 77, 1901 (1996).
  * (21) C. Micheletti, F. Seno, A. Maritan, and J. R. Banavar, Phys. Rev. Lett. 80, 2237 (1998).
  * (22) A. Irbäck, C. Peterson, F. Potthast, and E. Sandelin, Phys. Rev. E 58, R5249 (1998).
  * (23) A. Irbäck, C. Peterson, F. Potthast, and E. Sandelin, Structure 7, 347 (1999).
  * (24) Y. Iba, K. Tokita, and M. Kikuchi, J. Phys. Soc. Jpn. 67, 3985 (1998).
  * (25) K. Tokita, M. Kikuchi, and Y. Iba, Prog. Theor. Phys. Suppl. 138, 378 (2000).
  * (26) A. Rossi, A. Maritan, and C. Micheletti, J. Chem. Phys. 112, 2050 (2000).
  * (27) M. R. Betancourt and D. Thirumalai, J. Phys. Chem. B 106, 599 (2002).
  * (28) J. Zou and J. G. Saven. J. Chem. Phys. 118, 3843 (2003).
  * (29) Wang, B. Wang, Y. Liu, W. Chen, and C. Wang, Chin. Sci. Bull. 49, 426 (2004).
  * (30) X. Jiao, B. Wang, J. Su, W. Chen, and C. Wang, Phys. Rev. E 73, 061903-1 (2006).
  * (31) C. L Kleinman, N. Rodrigue, C. Bonnard, H. Philippe, and N. Lartillot, BMC Bioinform. 7, 326 (2006).
  * (32) I. Coluzza, H. G. Muller, and D. Frenkel, Phys. Rev. E 68, 046703-1 (2003).
  * (33) I. Coluzza, PLoS ONE, 6, e20853 (2011).
  * (34) A. W. Senior, R. Evans, J. Jumper, J. Kirkpatrick, L. Sifre, T. Green, C. Qin, A. Žídek, A. W. R. Nelson, A. Bridgland, H. Penedones, S. Petersen, K. Simonyan, S. Crossan, P. Kohli, D. T. Jones, D. Silver, K. Kavukcuoglu, and D. Hassabis, Nature 577, 706 (2020).
  * (35) Z. Li, Y. Yang, E. Faraggi, J. Zhan, and Y. Zhou, Proteins 82, 2565 (2014).
  * (36) J. Wang, H. Cao, J. Z. H. Zhang, and Y. Qi, Sci. Rep. 8, 6349 (2018).
  * (37) J. G. Greener, L. Moffat, and D. T. Jones, Sci. Rep. 8, 16189 (2018).
  * (38) D. Merk, L. Friedrich, F. Grisoni, and G. Schneider. Mol. Inf. 37, 1700153 (2018).
  * (39) A. Gupta, A. T. Müller, B. J. H. Huisman, J. A. Fuchs, P. Schneider, and G. Schneider, Mol. Inf. 37, 1700111 (2018).
  * (40) J. O’Connell, Z. Li, J. Hanson, R. Heffernan, J. Lyons, K. Paliwal, A. Dehzangi, Y. Yang, and Y. Zhou, Proteins 86, 629 (2018).
  * (41) A. Button, D. Merk, J. A. Hiss, and G. Schneider, Nat. Mach. Intell. 1, 307 (2019).
  * (42) F. C. Bernstein, T. F. Koetzle, G. J. B. Williams, E. F. Meyer Jr, M. D. Brice,J. R. Rodgers, O. Kennard, T. Shimanouchi, and M. Tasumi, Euro. J. Biochem. 80, 319 (1977).
  * (43) H. Li, R. Helling, C. Tang, and N. Wingreen, Science 273, 666 (1996).
  * (44) P. E. Wright and H. J. Dyson, J. Mol. Bio. 293, 321 (1999).
  * (45) V. N. Uversky, J. R. Gillespie, and A. L. Fink, Proteins 41, 415 (2000).
  * (46) A. K. Dunker, J. D. Lawson, C. J. Brown, R. M Williams, P. Romero, J. S. Oh, C. J. Oldfield, A. M. Campen, C. M. Ratliff, K. W. Hipps, J. Ausio, M. S. Nissen, R. Reeves, C. Kang, C. R. Kissinger, R. W. Bailey, M. D. Griswold, W. Chiu, E. C. Gerner, and Z. Obradovic, J. Mol. Graph. Model. 19, 26 (2001).
  * (47) P. Tompa, Trends Biochem. Sci. 27, 527 (2002).
  * (48) V. N. Uversky, Curr. Opin. Struct. Biol. 44, 18 (2017).