XXX XXX Original Paper

Wu et al.

[$\ast$]Corresponding author.

0Year 0Year 0Year

# Heterogeneous Entity Representation for Medicinal Synergy Prediction

Jiawei Wu  Jun Wen  Mingyuan Yan  Anqi Dong  Can Chen  School of Medicine,
National University of Singapore, Singapore 119077  Harvard Medical School,
Harvard University, Boston, MA 02115, USA  VA Boston Healthcare System, MA
02115, USA  Department of Mechanical and Aerospace Engineering, University of
California, Irvine, Irvine, CA 92697, USA  School of Data Science and Society,
University of North Carolina at Chapel Hill, Chapel Hill, NC 27599, USA
Department of Mathematics, University of North Carolina at Chapel Hill, Chapel
Hill, NC 27599, USA

(Date; Date; Date)

###### Abstract

Motivation: Medicinal synergy prediction is a powerful tool in drug discovery
and development that harnesses the principles of combination therapy to
enhance therapeutic outcomes by improving efficacy, reducing toxicity, and
preventing drug resistance. While a myriad of computational methods has
emerged for predicting synergistic drug combinations, a large portion of them
may overlook the intricate, yet critical relationships between various
entities in drug interaction networks, such as drugs, cell lines, and
diseases. These relationships are complex and multidimensional, requiring
sophisticated modeling to capture nuanced interplay that can significantly
influence therapeutic efficacy.
Results: We introduce a salient deep hypergraph learning method, namely,
Heterogeneous Entity Representation for MEdicinal Synergy prediction (HERMES),
to predict anti-cancer drug synergy. HERMES integrates heterogeneous data
sources, encompassing drug, cell line, and disease information, to provide a
comprehensive understanding of the interactions involved. By leveraging
advanced hypergraph neural networks with gated residual mechanisms, HERMES can
effectively learn complex relationships/interactions within the data. Our
results show HERMES demonstrates state-of-the-art performance, particularly in
forecasting new drug combinations, significantly surpassing previous methods.
This advancement underscores the potential of HERMES to facilitate more
effective and precise drug combination predictions, thereby enhancing the
development of novel therapeutic strategies.
Availability: The source code is publicly available at
https://github.com/Christina327/HERMES.
Contact<EMAIL_ADDRESS>

##  1 Introduction

Exploring drug combinations leads a promising avenue for enhancing cancer
treatment efficacy while minimizing toxicity and adverse reactions in modern
medicine (Jia et al.,, 2009; Csermely et al.,, 2013). Combination therapies,
involving multiple drugs administered as a single treatment regimen, offer
potential benefits over traditional single-drug approaches, particularly under
cancer and tumor treatment contexts (Chou,, 2006; O’Neil et al.,, 2016). Not
only do they hold the promise of greater therapeutic efficacy, but present an
opportunity to mitigate host toxicity and unwanted side effects, as the doses
of drug combinations are often sub-mutagenic compared to individual drug
doses. However, optimizing drug combinations can be challenging, as poorly
chosen combinations may lead to adverse effects and sub-optimal outcomes
(Hecht et al.,, 2009; Tol et al.,, 2009). Thus, there is a critical need to
identify precise synergistic drug pairs tailored to different cancer types.

Historically, the identification of effective combination drugs relied on
clinical experience, a process that is not only time-consuming and resource-
intensive but also prone to trial and error. In contrast, high-throughput
screening has emerged as an affordable and efficient strategy for identifying
synergistic drug pairs, leading to the generation of extensive datasets
(O’Neil et al.,, 2016; Holbeck et al.,, 2017; Jaaks et al.,, 2022). However,
certain limitations persist, such as the inability of cancer cell lines to
accurately represent in vivo states and the impracticality of exhaustively
testing all members of the full combinatorial space with high-throughput
screening (Ferreira et al.,, 2013; Goswami et al.,, 2015; Morris et al.,,
2016).

Recently, numerous computational methods for predicting drug synergy have been
proposed. Pioneering methods, including DeepSynergy (Preuer et al.,, 2018) and
Matchmaker (Kuru et al.,, 2021), utilize deep neural networks with both the
chemical properties of drugs and the gene expression profiles of cell lines.
The deep tensor factorization model (Sun et al.,, 2020) combines tensor
decomposition with neural networks to forecast the synergistic effects of drug
combinations. Additionally, TransSynergy (Liu and Xie,, 2021) adopts a
transformer network model using drug-target and protein-protein interaction
data. DeepDDS (Wang et al.,, 2022) employs graph convolutional networks (GCNs)
(Wu et al.,, 2020) and multi-layer perceptrons (MLP) for synergy prediction.
The current state-of-the-art method HypergraphSynergy (Liu et al.,, 2022) made
strides in this direction by incorporating hypergraph neural networks (HGNNs)
(Feng et al.,, 2019; Bai et al.,, 2021) to model these dynamics in a more
interconnected and multifaceted manner. Hypergraphs generalize graphs by
allowing hyperedges to connect more than two nodes (Chen and Rajapakse,, 2020;
Chen et al.,, 2021; Pickard et al.,, 2023), encode multidimensional (or high-
order) correlations and connections. However, HypergraphSynergy has notable
limitations in fully capturing the intricate nature of these interactions.
More importantly, its ability to generalize, especially in the context of
novel drug combinations, remains a challenge.

Figure 1: Overview of HERMES framework. HERMES has three key phases: (1)
initialization: acquiring and transforming initial features of drugs, cell
lines, and disease indications to a uniform dimensionality; (2) refinement:
enhancing feature representations through the construction of a dual-
relationship hypergraph and the application of hypergraph neural networks with
gated residual connections; (3) consolidation: integrating refined features
using a binary classification model to predict drug synergies with high
accuracy. Each phase is integral to the framework’s ability to process diverse
data types and generate precise drug synergy predictions.

In this article, we propose a novel deep hypergraph learning method –
Heterogeneous Entity Representation for MEdicinal Synergy prediction (HERMES)
– to enhance drug synergy prediction. HERMES distinguishes itself with its
innovative strategy for integrating a variety of data sources, including drug
chemical properties, cell line gene expressions, and interactions between
drugs and indications. This integration is achieved through a heterogeneous
hypergraph structure, which aids in assimilating extensive prior knowledge and
enhances the model’s ability to generalize in novel contexts. Another notable
breakthrough of HERMES is addressing the widespread issue of over-smoothing in
message-passing networks (Li et al.,, 2018) by implementing a gated residual
mechanism, which not only retains more information but also notably enhances
the expressiveness of the network. Our empirical results highlight HERMES’s
effectiveness, especially in novel scenarios, significantly surpassing
HypergraphSynergy and establishing it as a leading solution. The key
contributions of this article are (1) integration of varied knowledge sources
(including drugs, cell lines, and diseases) into a scalable and heterogeneous
model architecture; (2) enhancement of message passing over hypergraphs with
gated residual mechanisms for augmented network expressiveness; (3) empirical
evidence showcasing superior performance compared to previous methods, with a
particular emphasis on robust generalization in novel contexts.

This article is structured into four sections. The main architecture of HERMES
is detailed in Section 2. We assess the performance of HERMES (along with
other representative drug synergy prediction methods) including a
comprehensive ablation study in Section 3. Finally, we conclude by discussing
future research directions in Section 4.

##  2 Materials and Methods

###  2.1 Overview of HERMES

Our framework conceptualizes drug synergy prediction under a hypergraph
framework, treating it as a hyperedge prediction problem, where drugs, cell
lines, and diseases are represented by nodes while synergistic drug–drug–cell
line triplets and drug-disease pairwise relations are represented by
hyperedges. Hyperedge prediction is a generalization of edge prediction on
graphs (Chen and Liu,, 2023; Chen et al.,, 2023). This innovative approach
unfolds through three interconnected phases: initialization, refinement, and
consolidation (Figure 1). Each phase plays a pivotal role in processing and
integrating varied data types, thereby synthesizing comprehensive information
crucial for accurate drug synergy predictions. Drawing inspiration from
HypergraphSynergy, our model introduces novel methodologies and integrates
cutting-edge techniques to significantly enhance prediction accuracy and
reliability at each phase of analysis.

###  2.2 Feature Initialization

The initialization phase is crucial in our synergy prediction methodology, as
it involves acquisition of initial features for drugs, cell lines, and
diseases. This phase is pivotal in not only gathering diverse data but
ensuring these features are transformed to a uniform dimensionality, which is
fundamental for the subsequent analyses in our hypergraph model.

Drug features. We generate drug features by utilizing molecular graphs, which
are built from the Simplified Molecular Input Line Entry System (SMILES)
representations of drugs (Kim et al.,, 2019). By employing the
ConvMolFeaturizer method from the DeepChem library (Duvenaud et al.,, 2015),
we construct a molecular graph for each drug, with atoms represented as nodes
and bonds as edges. Subsequently, we employ advanced graph transformer
networks (GTNs) on molecular graphs to enhance the representations of atoms
and drugs. Denote the feature vector of atom $i$ as $\textbf{a}_{i}$. The
updated feature (denoted as $\textbf{a}^{\prime}_{i}$) through a GTN is then
computed as

$\displaystyle\textbf{a}^{\prime}_{i}=\sigma\Big{(}\textbf{W}_{1}\textbf{x}_{i}+\sum_{j\in\mathcal{N}(i)}\alpha_{ij}\textbf{W}_{2}\textbf{a}_{j}\Big{)},$
(1)

where $\mathcal{N}(i)$ denotes the neighboring atoms of atom $i$,
$\textbf{W}_{1}$ and $\textbf{W}_{2}$ are learnable weight matrices, and
$\sigma$ is a nonlinear activation function. Here, $\alpha_{ij}$ is the multi-
head attention coefficient defined as

$\alpha_{ij}=\text{Softmax}\Big{(}\frac{(\textbf{W}_{3}\textbf{a}_{i})^{\top}(\textbf{W}_{4}\textbf{a}_{j})}{\sqrt{d}}\Big{)},$

where $d$ is the latent size of each head, and $\textbf{W}_{3}$ and
$\textbf{W}_{4}$ are learnable weight matrices. GTNs can effectively capture
and model complex relationships and interactions within graph-structured data
by assigning attention scores to neighboring nodes, enabling a more targeted
aggregation of information compared to traditional convolutional methods such
as GCNs (Shi et al.,, 2020). Finally, we compute the feature vector of a drug
by aggregating all its updated atom features using maximum pooling. The choice
of maximum pooling is motivated by the fact that drug molecular graphs may
vary in the number of atom nodes. The resulting feature matrix for all drugs
is denoted as $\textbf{X}_{\text{drug}}$.

Cell line features. We utilize gene expression profiles to generate cell line
features. Following a log2 transformation and z-score normalization of gene
expression data, we employ an MLP to derive the features for cell lines,
ensuring uniform dimensionality with drug features. Denote the feature vector
for cell line $i$ by $\textbf{c}_{i}$, the updated feature (denoted as
$\textbf{c}^{\prime}_{i}$) through an MLP is computed as

$\textbf{c}^{\prime}_{i}=\sigma(\textbf{W}_{5}\textbf{c}_{i}+\textbf{b}_{5}),$
(2)

where $\textbf{W}_{5}$ and $\textbf{b}_{5}$ are learnable weight matrix and a
bias vector, respectively. The resulting feature matrix for all cell lines is
denoted as $\textbf{X}_{\text{cell}}$.

Disease features. We create disease features by leveraging a novel embedding
method called CODER, a knowledge-graph-guided large language model for cross-
lingual medical term representation using contrastive learning (Yuan et al.,,
2022). CODER excels in providing close vector representations for various
terms representing similar medical concepts in multiple languages, making it
particularly beneficial for extracting enriched and contextually relevant
disease features (Yuan et al.,, 2022). Similar to cell line features
processing, we employ an MLP on the obtained disease features to align with
the consistent dimensionality of drug and cell line features. The resulting
feature matrix for all diseases is denoted as $\textbf{X}_{\text{dis}}$.
Notably, integrating CODER with disease information enhances our model’s
ability to comprehend the intricate relationships between drugs, cell lines,
and diseases within the hypergraph framework.

###  2.3 Feature Refinement

In the refinement phase, we focus on transforming the obtained drug, cell
line, and disease features into contextually enriched representations. This
process is achieved through the construction of a dual-relation hypergraph and
the implementation of an advanced feature enhancement technique, namely HGNNs
with gated residual connections. These elements collectively elevate our
model’s ability to identify complex interrelationships, setting the stage for
more precise synergy predictions.

Hypergraph construction. We construct a novel dual-relationship hypergraph
involving relationships between drugs, cell lines, and diseases. The
hypergraph contains two types of hyperedges: (1) drug-drug-cell line triplets
(third-order interactions), which capture the potential synergistic effects
between specific drugs and cell lines; (2) drug-disease relationship (pairwise
interactions or traditional edges), which represents associations between
drugs and their indications (medical conditions they treat), providing
insights into the therapeutic spectrum of drugs. This dual-relationship
hypergraph, which has never been considered in previous methods, is pivotal in
enriching the feature representations by integrating the various features of
drugs, cell lines, and diseases.

Hypergraph neural networks. We employ advanced HGNNs with gated residual
mechanisms to refine drug, cell lines, and disease features. HGNNs excel in
capturing the complex relationships and interactions within the hypergraph,
making it a pivotal component in our model for drug synergy prediction (Feng
et al.,, 2019). Denote the feature matrix as
$\textbf{X}=\text{vercat}(\textbf{X}_{\text{drug}},\textbf{X}_{\text{cell}},\textbf{X}_{\text{dis}})$
where vercat is the vertical concatenation operation. The refined features
(denoted as $\textbf{X}^{\prime}$) through an HGNN are computed as

$\textbf{X}^{\prime}=\sigma(\textbf{D}^{-1}\textbf{H}\textbf{E}^{-1}\textbf{H}^{\top}\textbf{X}\textbf{W}_{6}),$
(3)

where D and E are diagonal matrices of node and hyperedge degrees,
respectively. H is the incidence matrix of the hypergraph, containing binary
values to indicate the presence or absence of any node in any hyperedge, and
$\textbf{W}_{6}$ is learnable weight matrix. Here, the superscripts ${-1}$ and
$\top$ denote matrix inversion and transposition, respectively. To enhance the
expressive capacity of our network and extract more intricate patterns from
the hypergraph, we increase the number of layers in the HGNN. However, this
approach leads to a well-known challenge in graph neural networks, i.e., over-
smoothing (Cai and Wang,, 2020).

The phenomenon of over-smoothing occurs when the network layers become too
deep, leading to the homogenization of node features and a loss of valuable
information. To mitigate this issue and further empower our network’s
expressive capabilities, we introduce a novel solution using gated residual
connections (Li et al.,, 2018). The implementation of gated residual
connections within an HGNN is defined as

$\textbf{X}^{\prime}=\textbf{X}+\sigma\Big{(}\textbf{W}_{7}\sigma(\textbf{D}^{-1}\textbf{H}\textbf{E}^{-1}\textbf{H}^{\top}\textbf{X}\textbf{W}_{6})+\textbf{b}_{7}\Big{)}\textbf{X},$
(4)

where $\textbf{W}_{7}$ and $\textbf{b}_{7}$ are learnable weight matrix and
bias vector of the gate, respectively. The introduction of the gated residual
connection in HGNNs allows our neural network to dynamically regulate the
integration of original and convoluted features, which significantly reduces
the risk of over-smoothing by providing a controlled blending of features.
Additionally, we implement the equilibrium bias initialization (EBI) strategy
(Wang et al.,, 2010), where we initialize $\textbf{b}_{7}$ with a relatively
negative value. The EBI strategy ensures that the gate functions close to an
identity operation at the start of training, allowing for gradual and more
effective feature integration. Our refinement approach utilizing HGNNs coupled
with gated residual connections significantly enhances the model’s ability to
discern complex patterns and interactions between drugs, cell lines, and
diseases, thereby creating more accurate representations for them compared to
previous methods such as HypergraphSynergy.

###  2.4 Feature Consolidation

In the consolidation phase, our primary objective is to produce predictive
insights by integrating the refined features obtained from the earlier phase.
We employ a binary classification model that uses the features refined over
the dual-relationship hypergraph to identify drug synergies. It classifies the
interactions between drug combinations and cell lines into two categories:
synergistic and non-synergistic, which provides a clear, binary output for
each potential drug synergy scenario.

Predictive model. The essence of the consolidation phase is using a predictive
model that leverages the features obtained in the refinement phase to evaluate
the potential synergistic effects of drug combinations on specific cell lines
and their implications for disease treatment. Suppose the final refined
features of drug $i$, drug $j$, and cell line $k$ are denoted by
$\tilde{\textbf{d}}_{i}$, $\tilde{\textbf{d}}_{j}$, and
$\tilde{\textbf{c}}_{k}$, respectively. The synergy score (denoted as $S$) can
be computed using an MLP, i.e.,

$S=\sigma\Big{(}\textbf{W}_{8}\text{vercat}(\tilde{\textbf{d}}_{i},\tilde{\textbf{d}}_{j},\tilde{\textbf{c}}_{k})+\textbf{b}_{8}\Big{)},$
(5)

where $\textbf{W}_{8}$ and $\textbf{b}_{8}$ are a learnable weight matrix and
a bias vector, respectively.

To address the unordered nature of drug pairs, a critical aspect in synergy
prediction, we implement a pairwise symmetric permutation augmentation
strategy (Zhou et al.,, 2024). This approach involves presenting both
permutations of each drug pair in our dataset. For example, if a data sample
includes the combination (drug1, drug2, cell line), we also introduce (drug2,
drug1, cell line) as a distinct sample. This augmentation is essential to
ensure that our model is agnostic to the order of drugs, enhancing its
capability to uniformly recognize synergy.

Model training. Given the binary nature of our prediction task, we utilize the
cross-entropy loss function (Mao et al.,, 2023), a standard and effective
choice for binary classification models defined as

$\text{Loss}=-\frac{1}{N}\sum_{i=1}^{N}[y_{i}\log(\hat{y}_{i})+(1-y_{i})\log(1-\hat{y}_{i})],$
(6)

where $y_{i}$ and $\hat{y}_{i}$ represent the actual and predicted label
(synergistic or non-synergistic) of the $i$th sample, respectively, and $N$ is
the total number of samples. The model is trained on a dataset consisting of
drug pairs, cell lines, and their corresponding synergy labels. The training
process involves adjusting the weights of the neural network to minimize the
cross-entropy loss, enhancing the model’s ability to accurately classify drug
synergies.

###  2.5 HERMES and HypergraphSynergy

Although HERMES is inspired by HypergraphSynergy, it significantly diverges
from HypergraphSynergy in the following aspects:

  * •

HERMES integrates an additional knowledge source of diseases and constructs a
dual-relationship hypergraph.

  * •

HERMES leverages advanced GTNs to obtain drug features, while
HypergraphSynergy only utilizes traditional GCNs.

  * •

HERMES improves HGNNs during the refinement phase by incorporating gated
residual mechanisms to address the issue of over-smoothing.

  * •

HERMES exploits useful learning techniques such as EBI and pairwise symmetric
permutation augmentation to enhance the model’s learning ability.

Below, we show that HERMES significantly outperforms HypergraphSynergy in two
drug synergy datasets, especially under the context of predicting new drug
combinations.

##  3 Experiments

###  3.1 Datasets

We collected four categories of data, including drug synergy data, molecular
information for drugs, genomic characteristics of cancer cell lines, and
indications for drug therapy in diseases, from multiple publicly available
databases. The details of these data are listed as follows:

  * •

Drug synergy datasets: We gathered data on the synergy of anti-cancer drugs
from two prominent large-scale tumor screening datasets - the O’Neil dataset
(O’Neil et al.,, 2016) and the NCI-ALMANAC dataset (Holbeck et al.,, 2017).
The O’Neil dataset comprises $23,062$ samples involving 38 unique drugs and
$39$ distinct human cancer cell lines. Each sample measures Loewe synergy
scores for two drugs in combination with a specific cell line. The NCI-ALMANAC
dataset contains $304,549$ samples, including ComboScores for $104$ FDA-
approved drugs in pairings across the NCI-60 cell line panel.

  * •

Drug molecular structures: Information on the SMILES of drugs is obtained from
the PubChem database (Kim et al.,, 2019).

  * •

Gene expression in cancer cell lines: Data on gene expression in cancer cell
lines are sourced from the Cell Lines Project within the COSMIC database
(Forbes et al.,, 2015). In this context, we specifically considered the
expression values of $651$ genes related to the COSMIC cancer gene census.
These expression values are subjected to logarithmic (log2) transformation and
z-score normalization.

  * •

Indications for drug therapy in diseases: The indication dataset is sourced
from PrimeKG (Chandak et al.,, 2023), a multi-modal knowledge graph designed
to support precision medicine. PrimeKG integrates data from $20$ high-quality
biomedical resources, including DisGeNET and DrugBank. It provides detailed
relationships between drugs and diseases, covering indications,
contraindications, and off-label uses. This dataset captures how drugs target
disease-associated molecular perturbations, facilitating analyses in drug
repurposing and therapeutic strategies

We performed a comprehensive data preprocessing on two primary drug synergy
datasets – NCI-ALMANAC and O’Neil. To ensure the quality and relevance of the
data, we excluded cell lines lacking gene expression information and drugs
without SMILES details. Following this preprocessing phase, the NCI-ALMANAC
dataset comprises $74,139$ measurement samples of ComboScores for $87$ drugs
across $55$ cancer cell lines, while the O’Neil dataset encompasses $18,950$
samples of Loewe synergy scores for $38$ drugs and $39$ cancer cell lines
(Table 1). Subsequently, we removed drugs from the ‘drug-disease dataset’ that
are not present in the aforementioned drug synergy datasets. This led to the
extraction of indications for $111$ diseases corresponding to $111$ drugs
within the NCI-ALMANAC dataset and $222$ diseases corresponding to $222$ drugs
within the O’Neil dataset.

Table 1: Statistics of the two datasets. Dataset | #Drugs | #Targets | #Samples
---|---|---|---
NCI-ALMANAC | 87 | 55 | 74139
O’Neil | 38 | 39 | 18950

###  3.2 Baselines

We conducted a comparative analysis of our approach with representative drug
synergy prediction techniques. Below is a brief overview of each of the
baseline methods:

  * •

DeepSynergy (Preuer et al.,, 2018): It utilizes a three-layer feedforward
neural network to predict synergy scores, incorporating gene expression as
cell line features and three types of chemical descriptors as drug features.

  * •

DTF (Sun et al.,, 2020): It extracts latent features from the drug synergy
matrix through tensor factorization and employs them to train a deep neural
network model for predicting drug synergy.

  * •

HypergraphSynergy (Liu et al.,, 2022) (the current state-of-the-art method):
It formulates synergistic drug combinations across cancer cell lines as a
hypergraph. In this hypergraph, drugs and cell lines are represented by nodes,
while synergistic drug–drug–cell line triplets are represented by hyperedges.
It leverages the biochemical features of drugs and cell lines as node
attributes. Additionally, a HGNN is employed to learn drug and cell line
embeddings from the hypergraph and predict drug synergy.

###  3.3 Experiment Setup

In this study, each dataset is initially divided into two distinct sets: a
training set, accounting for $90\%$ of the total data, and a test set
comprising the remaining $10\%$. The training set undergoes a rigorous five-
fold cross-validation process. This process is structured in three unique
partitioning strategies to ensure comprehensive evaluation:

  * •

Random: Samples are randomly divided, providing a baseline assessment of model
performance.

  * •

Target: Samples are stratified by target cell line, ensuring each fold’s
validation set contains unique cell lines not present in its training set.

  * •

DrugComb: Samples are stratified based on drug combinations. Each validation
set included drug combinations not seen in the training set, although
individual drugs might overlap.

The test set is used for final evaluation, whereby an unbiased assessment of
the model’s predictive power is ensured.

Additionally, for the classification task, synergy scores are converted to
binary outcomes. Following established protocols (Preuer et al.,, 2018; Sun et
al.,, 2020), a threshold of 30 is used. Scores above this threshold indicate a
positive synergy, while scores below are deemed negative. Ultimately, the
model’s performance is evaluated using key metrics: area under the receiver
operating characteristic curve (AUROC) and area under the precision-recall
curve (AUPRC). Both metrics provide comprehensive view of the model’s
effectiveness in precision and recall aspects.

###  3.4 Hyperparameters

To optimize the performance of our model, we conducted a grid search to fine-
tune hyperparameters such as the refinement layer and learning rate. Table 2
presents the range of values considered for each hyperparameter. The optimal
values (in bold) are selected based on their performance during the training
and validation phases.

Table 2: Hyperparameter selection (Selected hyperparameters are highlighted in bold.) Hyperparameter | Values
---|---
learning rate | {1e-3, 5e-4, 2e-4, 1e-4, 5e-5, 2e-5}
weight decay | {1e-1, 1e-2, 1e-3, 1e-4}
attention heads | {2, 4, 8}
refinement layer | {2, 3, 4}
interaction weight | {0, 0.02, 0.05, 0.1, 0.2, 0.5, 1.0}

###  3.5 Performance Comparison and Analysis

RandomTargetDrugComb$65$$70$$75$$80$$85$$90$***AUROC(%)HERMESHypergraphSynergyDTFDeepSynergyARandomTargetDrugComb$25$$30$$35$$40$$45$$50$$55$$60$*AUPRC(%)

RandomTargetDrugComb$80$$82$$84$$86$$88$$90$$92$$94$$96$****AUROC(%)BRandomTargetDrugComb$30$$40$$50$$60$$70$***AUPRC(%)

Figure 2: Model performance comparison for (A) ALMANAC Dataset and (B) ONeil
Dataset. The left panels shows the AUROC (%) scores for different models
across three validation modes (Random, Target, DrugComb). The right panels
present the AUPRC (%) scores for the same models and validation modes. (B)
ONeil Dataset. Red asterisks indicate statistical significance between HERMES
and HypergraphSynergy(*** $p$-value $<0.001$; ** $p$-value $<0.01$; *
$p$-value $<0.05$; two-sample $t$-test).

In the NCI-ALMANAC dataset, HERMES exhibits superior performance across all
evaluation strategies compared to the baseline methods (Figure 2A). Under the
random partitioning strategy, HERMES achieves an AUROC of $85.91\%$, which is
significantly higher than HypergraphSynergy, DTF, and DeepSynergy, which score
$85.30\%$, $82.38\%$, and $83.50\%$, respectively (HERMES vs.
HypergraphSynergy: $p$-value $<0.05$, two-sample $t$-test). In the more
challenging ‘DrugComb’ mode, HERMES achieves an AUROC of $79.75\%$,
significantly outperforming the other methods, with HypergraphSynergy
recording $77.98\%$ ($p$-value $<0.01$, two-sample $t$-test), and DTF and
DeepSynergy achieving even lower scores. For the ‘Target’ mode, HERMES still
performs better than the rest of the methods. Additionally, the AUPRC results
reflect a similar trend (right panel).

The performance of HERMES is further validated using the ONEIL dataset (Figure
2B), and similar results are observed across the three modes in terms of AUROC
and AUPRC. In the random stratification, HERMES attains an AUROC of 93.67%,
significantly outperforming HypergraphSynergy, DTF, and DeepSynergy, which
score 92.30%, 91.38%, and 90.60%, respectively (HERMES vs. HypergraphSynergy:
$p$-value $<0.001$, two-sample $t$-test). In the ‘DrugComb’ mode, HERMES
maintains a competitive edge with an AUROC of 88.34%, significantly surpassing
HypergraphSynergy at 86.22% ($p$-value $<0.05$, two-sample $t$-test). HERMES
also achieves better performance under the ‘Target’ mode compared to the other
methods.

DrugSingleDrugDouble$55$$60$$65$$70$$75$******AUROC(%)DrugSingleDrugDouble$16$$18$$20$$22$$24$$26$$28$$30$**AUPRC(%)
Figure 3: Performance comparison of HERMES and HypergraphSynergy in the
DrugSingle and DrugDouble modes using the ALMANAC dataset. Red asterisks
indicate statistical significance between HERMES and HypergraphSynergy(***
$p$-value $<0.001$; ** $p$-value $<0.01$; two-sample $t$-test).

To further validate the performance of our model, we considered two new
validation modes using the ALMANAC dataset (which contains more samples) to
test the model’s generalization ability to previously unseen drugs. The two
modes are defined as follows: (1) DrugSingle where samples are stratified
based on individual drugs, with one drug in each drug combination in the
validation set being novel to the training set. This approach reduces the
dataset size and challenges the model to predict synergies involving new
drugs; (2) DrugDouble where samples are stratified such that each drug
combination in the validation set contains two drugs, neither of which appears
in the training set, further reducing the dataset size and increasing the
difficulty. The results of ‘DrugSignle’ and ‘DrugDouble’ are shown in Figure
3. In the ‘DrugSingle’ mode, HERMES achieves an AUROC of 72.13%, which is
extremely significantly higher than HypergraphSynergy’s 67.84% ($p$-value
$<0.001$, two-sample $t$-test). Similarly, in the ‘DrugDouble’ mode, HERMES
scores an AUROC of 68.39%, highly significantly outperforming
HypergraphSynergy ($p$-value $<0.001$, two-sample $t$-test). The AUPRC metrics
also show a consistent advantage for HERMES.

The comprehensive evaluation across different datasets and validation modes
underscores the efficacy of HERMES in drug synergy prediction. The model
consistently outperforms baseline methods in both AUROC and AUPRC metrics,
demonstrating robust performance even in challenging scenarios. The results
validate the superiority of HERMES in handling diverse drug synergy prediction
tasks, emphasizing its importance for predicting untested drug synergies in
clinical applications.

###  3.6 Ablation Study

The ablation study, as detailed in Table 3, scrutinizes the contributions of
various components within the HERMES model utilizing the NCI-ALMANAC dataset.
The evaluation encompasses three distinct modes: ‘Random’, ‘Target’, and
‘DrugComb’, with the primary metric being the average test AUROC for each
configuration. The benchmark performance of the complete HERMES model
registers an AUROC of $85.9\%$ in the ‘Random’ mode, $79.4\%$ in the ‘Target’
mode, and $79.8\%$ in the ‘DrugComb’ mode. These figures establish the
baseline for assessing the impact of systematically removing specific
components.

Table 3: Ablation study results. Average test AUROC of HERMES with and without key components. Strategy | Random | Target | DrugComb
---|---|---|---
HERMES | 85.9% | 79.4% | 79.8%
w/o transformer | 82.8% | 78.5% | 78.3%
w/o disease | 85.7% | 79.2% | 79.5%
w/o gated residual | 84.0% | 79.0% | 78.8%
w/o gate | 82.5% | 78.7% | 78.5%
w/o weighted hyperedge | 82.8% | 78.6% | 78.5%

Contribution of transformer architecture. Excluding the transformer component
results in significant performance drops across all three modes. These
findings highlight the indispensable role of the transformer in capturing
intricate relationships within the data.

Impact of disease knowledge. The exclusion of disease knowledge demonstrates a
slight yet notable impact on the model’s performance. These findings suggest
that incorporating disease knowledge provides a performance advantage, but the
margin is relatively narrow. Future work could explore more sophisticated
methods of integrating disease knowledge to potentially amplify this
performance benefit.

Role of gated residuals. Omitting the gated residual connections results in a
notable performance degradation, particularly in the ‘Random’ mode. This
emphasizes the significance of gated residuals in maintaining high predictive
accuracy, especially in scenarios with higher data variability.

Importance of gating mechanisms. The removal of the gating mechanism induces
the most pronounced decline in performance across all modes. This underscores
the critical function of gating mechanisms in modulating information flow and
enhancing the model’s resilience.

Effect of weighted hyperedges. The absence of weighted hyperedges also leads
to considerable performance decline across all three modes. This underscores
the necessity of weighted hyperedges for accurately modeling complex
interactions, particularly in the diverse Random mode.

In summary, the ablation study clearly demonstrates that each component of the
HERMES model contributes uniquely to its overall efficacy. Transformer
architecture, disease knowledge, gated residuals, gating mechanisms, and
weighted hyperedges collectively enhance the model’s predictive performance.
These results affirm the sophisticated design of HERMES, showcasing the
synergistic effect of its components in achieving superior predictive accuracy
across various validation strategies.

##  4 Conclusion

In this article, we introduced HERMES, a novel deep hypergraph learning method
aimed at enhancing drug synergy prediction. By incorporating disease
information and employing advanced GTNs and HGNNs with gated residual
connections, HERMES achieves a superior performance on two drug synergy
datasets compared to HypergraphSynergy and other methods. For future work, we
will consider improving the scope of our analysis by incorporating a wider
array of heterogeneous data sources such as target proteins, biological
pathways, and therapeutic markers. We aim to bridge the gap between
theoretical drug synergy prediction and practical clinical applications,
ultimately contributing to more effective and personalized treatment
strategies.

## References

  * Bai et al., (2021) Bai, S., Zhang, F., and Torr, P. H. (2021).  Hypergraph convolution and hypergraph attention.  Pattern Recognition, 110:107637.
  * Cai and Wang, (2020) Cai, C. and Wang, Y. (2020).  A note on over-smoothing for graph neural networks.  arXiv preprint arXiv:2006.13318.
  * Chandak et al., (2023) Chandak, P., Huang, K., and Zitnik, M. (2023).  Building a knowledge graph to enable precision medicine.  Scientific Data, 10(67).
  * Chen et al., (2023) Chen, C., Liao, C., and Liu, Y.-Y. (2023).  Teasing out missing reactions in genome-scale metabolic networks through hypergraph learning.  Nature Communications, 14(1):2375.
  * Chen and Liu, (2023) Chen, C. and Liu, Y.-Y. (2023).  A survey on hyperlink prediction.  IEEE Transactions on Neural Networks and Learning Systems.
  * Chen and Rajapakse, (2020) Chen, C. and Rajapakse, I. (2020).  Tensor entropy for uniform hypergraphs.  IEEE Transactions on Network Science and Engineering, 7(4):2889–2900.
  * Chen et al., (2021) Chen, C., Surana, A., Bloch, A. M., and Rajapakse, I. (2021).  Controllability of hypergraphs.  IEEE Transactions on Network Science and Engineering, 8(2):1646–1657.
  * Chou, (2006) Chou, T.-C. (2006).  Theoretical basis, experimental design, and computerized simulation of synergism and antagonism in drug combination studies.  Pharmacological reviews, 58(3):621–681.
  * Csermely et al., (2013) Csermely, P., Korcsmáros, T., Kiss, H. J., London, G., and Nussinov, R. (2013).  Structure and dynamics of molecular networks: a novel paradigm of drug discovery: a comprehensive review.  Pharmacology & therapeutics, 138(3):333–408.
  * Duvenaud et al., (2015) Duvenaud, D. K., Maclaurin, D., Iparraguirre, J., Bombarell, R., Hirzel, T., Aspuru-Guzik, A., and Adams, R. P. (2015).  Convolutional networks on graphs for learning molecular fingerprints.  Advances in neural information processing systems, 28.
  * Feng et al., (2019) Feng, Y., You, H., Zhang, Z., Ji, R., and Gao, Y. (2019).  Hypergraph neural networks.  In Proceedings of the AAAI conference on artificial intelligence, volume 33, No.01, pages 3558–3565.
  * Ferreira et al., (2013) Ferreira, D., Adega, F., and Chaves, R. (2013).  The importance of cancer cell lines as in vitro models in cancer methylome analysis and anticancer drugs testing.  Oncogenomics and cancer proteomics-novel approaches in biomarkers discovery and therapeutic targets in cancer, 1:139–166.
  * Forbes et al., (2015) Forbes, S. A., Beare, D., Gunasekaran, P., Leung, K., Bindal, N., Boutselakis, H., Ding, M., Bamford, S., Cole, C., Ward, S., et al. (2015).  Cosmic: exploring the world’s knowledge of somatic mutations in human cancer.  Nucleic acids research, 43(D1):D805–D811.
  * Goswami et al., (2015) Goswami, C. P., Cheng, L., Alexander, P., Singal, A., and Li, L. (2015).  A new drug combinatory effect prediction algorithm on the cancer cell based on gene expression and dose–response curve.  CPT: Pharmacometrics & Systems Pharmacology, 4(2):80–90.
  * Hecht et al., (2009) Hecht, J. R., Mitchell, E., Chidiac, T., Scroggin, C., Hagenstad, C., Spigel, D., Marshall, J., Cohn, A., McCollum, D., Stella, P., et al. (2009).  A randomized phase iiib trial of chemotherapy, bevacizumab, and panitumumab compared with chemotherapy and bevacizumab alone for metastatic colorectal cancer.  Journal of Clinical Oncology, 27(5):672–680.
  * Holbeck et al., (2017) Holbeck, S. L., Camalier, R., Crowell, J. A., Govindharajulu, J. P., Hollingshead, M., Anderson, L. W., Polley, E., Rubinstein, L., Srivastava, A., Wilsker, D., et al. (2017).  The national cancer institute almanac: a comprehensive screening resource for the detection of anticancer drug pairs with enhanced therapeutic activity.  Cancer research, 77(13):3564–3576.
  * Jaaks et al., (2022) Jaaks, P., Coker, E. A., Vis, D. J., Edwards, O., Carpenter, E. F., Leto, S. M., Dwane, L., Sassi, F., Lightfoot, H., Barthorpe, S., et al. (2022).  Effective drug combinations in breast, colon and pancreatic cancer cells.  Nature, 603(7899):166–173.
  * Jia et al., (2009) Jia, J., Zhu, F., Ma, X., Cao, Z. W., Li, Y. X., and Chen, Y. Z. (2009).  Mechanisms of drug combinations: interaction and network perspectives.  Nature reviews Drug discovery, 8(2):111–128.
  * Kim et al., (2019) Kim, S., Chen, J., Cheng, T., Gindulyte, A., He, J., He, S., Li, Q., Shoemaker, B. A., Thiessen, P. A., Yu, B., et al. (2019).  Pubchem 2019 update: improved access to chemical data.  Nucleic acids research, 47(D1):D1102–D1109.
  * Kuru et al., (2021) Kuru, H. I., Tastan, O., and Cicek, A. E. (2021).  Matchmaker: a deep learning framework for drug synergy prediction.  IEEE/ACM transactions on computational biology and bioinformatics, 19(4):2334–2344.
  * Li et al., (2018) Li, Q., Han, Z., and Wu, X.-M. (2018).  Deeper insights into graph convolutional networks for semi-supervised learning.  In Proceedings of the AAAI conference on artificial intelligence, volume 32, No.1.
  * Liu and Xie, (2021) Liu, Q. and Xie, L. (2021).  Transynergy: Mechanism-driven interpretable deep neural network for the synergistic prediction and pathway deconvolution of drug combinations.  PLoS computational biology, 17(2):e1008653.
  * Liu et al., (2022) Liu, X., Song, C., Liu, S., Li, M., Zhou, X., and Zhang, W. (2022).  Multi-way relation-enhanced hypergraph representation learning for anti-cancer drug synergy prediction.  Bioinformatics, 38(20):4782–4789.
  * Mao et al., (2023) Mao, A., Mohri, M., and Zhong, Y. (2023).  Cross-entropy loss functions: Theoretical analysis and applications.  arXiv preprint: 2304.07288.
  * Morris et al., (2016) Morris, M. K., Clarke, D. C., Osimiri, L. C., and Lauffenburger, D. A. (2016).  Systematic analysis of quantitative logic model ensembles predicts drug combination effects on cell signaling networks.  CPT: pharmacometrics & systems pharmacology, 5(10):544–553.
  * O’Neil et al., (2016) O’Neil, J., Benita, Y., Feldman, I., Chenard, M., Roberts, B., Liu, Y., Li, J., Kral, A., Lejnine, S., Loboda, A., et al. (2016).  An unbiased oncology compound screen to identify novel combination strategies.  Molecular cancer therapeutics, 15(6):1155–1162.
  * Pickard et al., (2023) Pickard, J., Chen, C., Salman, R., Stansbury, C., Kim, S., Surana, A., Bloch, A., and Rajapakse, I. (2023).  Hat: Hypergraph analysis toolbox.  PLOS Computational Biology, 19(6):e1011190.
  * Preuer et al., (2018) Preuer, K., Lewis, R. P., Hochreiter, S., Bender, A., Bulusu, K. C., and Klambauer, G. (2018).  Deepsynergy: predicting anti-cancer drug synergy with deep learning.  Bioinformatics, 34(9):1538–1546.
  * Shi et al., (2020) Shi, Y., Huang, Z., Feng, S., Zhong, H., Wang, W., and Sun, Y. (2020).  Masked label prediction: Unified message passing model for semi-supervised classification.  arXiv preprint: 2009.03509.
  * Sun et al., (2020) Sun, Z., Huang, S., Jiang, P., and Hu, P. (2020).  Dtf: deep tensor factorization for predicting anticancer drug synergy.  Bioinformatics, 36(16):4483–4489.
  * Tol et al., (2009) Tol, J., Koopman, M., Cats, A., Rodenburg, C. J., Creemers, G. J., Schrama, J. G., Erdkamp, F. L., Vos, A. H., van Groeningen, C. J., Sinnige, H. A., et al. (2009).  Chemotherapy, bevacizumab, and cetuximab in metastatic colorectal cancer.  New England Journal of Medicine, 360(6):563–572.
  * Wang et al., (2010) Wang, H.-H., Klessen, R. S., Dullemond, C. P., Van Den Bosch, F. C., and Fuchs, B. (2010).  Equilibrium initialization and stability of three-dimensional gas discs.  Monthly Notices of the Royal Astronomical Society, 407(2):705–720.
  * Wang et al., (2022) Wang, J., Liu, X., Shen, S., Deng, L., and Liu, H. (2022).  Deepdds: deep graph neural network with attention mechanism to predict synergistic drug combinations.  Briefings in Bioinformatics, 23(1):bbab390.
  * Wu et al., (2020) Wu, Z., Pan, S., Chen, F., Long, G., Zhang, C., and Philip, S. Y. (2020).  A comprehensive survey on graph neural networks.  IEEE transactions on neural networks and learning systems, 32(1):4–24.
  * Yuan et al., (2022) Yuan, Z., Zhao, Z., Sun, H., Li, J., Wang, F., and Yu, S. (2022).  Coder: Knowledge-infused cross-lingual medical term embedding for term normalization.  Journal of biomedical informatics, 126:103983.
  * Zhou et al., (2024) Zhou, A., Yang, K., Burns, K., Cardace, A., Jiang, Y., Sokota, S., Kolter, J. Z., and Finn, C. (2024).  Permutation equivariant neural functionals.  Advances in Neural Information Processing Systems, 36.


High quality factor mechanical resonance in a silicon nanowire High quality
factor mechanical resonance in a silicon nanowire Presnov, Kafanov, Dorofeev,
Bozhev, Trifonov, Pashkin, Krupenin Presnov, Kafanov, Dorofeev, Bozhev,
Trifonov, Pashkin, Krupenin *

# High quality factor mechanical resonance in a silicon nanowire

D. E. Presnov†,∘ S. Kafanov‡ A. A. Dorofeev† I. V. Bozhev† A. S. Trifonov†,∘
Yu. A<EMAIL_ADDRESS>V. A.
<EMAIL_ADDRESS>[[ †Quantum Technology Centre, Faculty of
Physics, M. V. Lomonosov Moscow State University, Moscow, 119991, Russia.
∘D. V. Skobeltsyn Institute of Nuclear Physics, M. V. Lomonosov Moscow State
University, Moscow, 119991, Russia.
‡Department of Physics, Lancaster University, Lancaster, LA1 4YB, United
Kingdom.
⋆Lebedev Physical Institute, Moscow 119991, Russia.

###### Abstract

Resonance properties of nanomechanical resonators based on doubly clamped
silicon nanowires, fabricated from silicon-on-insulator and coated with a thin
layer of aluminum, were experimentally investigated. Resonance frequencies of
the fundamental mode were measured at a temperature of $20\,\mathrm{mK}$ for
nanowires of various sizes using the magnetomotive scheme. The measured values
of the resonance frequency agree with the estimates obtained from the Euler-
Bernoulli theory. The measured internal quality factor of the $5\,\mathrm{\mu
m}$-long resonator, $3.62\times 10^{4}$, exceeds the corresponding values of
similar resonators investigated at higher temperatures. The structures
presented can be used as mass sensors with an expected sensitivity $\sim
6\times 10^{-20}\,\mathrm{g}\,\mathrm{Hz}^{-1/2}$.

62.25.-g, 81.07.-b, 81.07.Oj, 07.07.Df

## Introduction

Nanoelectromechanical systems have been widely used in the fundamental and
applied physics research [1, 2, 3, 4], due to the rapid development of
technology in recent years. Various microelectromechanical systems have become
an integral part of our daily lives. Such structures are present almost
everywhere, be it a phone or a car, in the form of various micron-scale
accelerometers, gyroscopes, etc. Today, when electronic devices make use of
the properties of individual atoms and molecules [5, 6, 7], applications of
nanoelectromechanical systems include ultrasensitive detection of mass, down
to the mass of single molecules [8, 9], force [10], pressure [11] and
displacement [12, 13]. By coupling nanomechanical resonators with optical and
electronic transducers, it became possible to explore various quantum effects
[14, 15, 16]. Also, small-size mechanical resonators are effective tools in
understanding properties of superfluids [17].

The operation principle of nanoelectromechanical mass sensors is based on the
phenomenon of mechanical resonance that occurs when oscillations are excited,
for example, in beams clamped at both ends, with a sub-micron cross-sectional
area. The eigenfrequency of the resonator depends on its geometry and
material, and increases when the resonator dimensions get smaller. It was
shown [18] that the quality factor of resonators depends on their surface-to-
volume ratio, and is also affected by the surface quality of the fabricated
structure. The change of the resonator mass leads to the shift of its
resonance frequency. The minimum detectable change of the mass is given by the
expression [19]:

$\delta m=-2\frac{0.735m}{f_{0}}\delta f_{0},$ (1)

where $m$ is the initial nanowire mass, $\delta m$ is the change of the mass,
$f_{0}$ is the resonance frequency, and $\delta f_{0}$ is the minimum
detectable frequency shift determined by the internal noise of the resonator
and noise of the measurement system. Thus, to build an ultra-sensitive mass
sensor, it is desirable to reduce the effective mass of the resonator and
increase its resonance frequency, while simultaneously increasing its quality
factor in order to resolve as small frequency shift as possible. Up to now,
the record value reported for the quality factor of silicon resonators with a
sub-micron cross-section was $\sim 1.8\times 10^{4}$ [20].

In the past decade, silicon-on-insulator (SOI) became a commonly used material
for the fabrication of silicon nanoelectromechanical systems [21]. It comes as
a three-layer wafer, in which the upper thin layer of single-crystal silicon
is separated from the base substrate by a thin inter-layer of the silicon
oxide. This material is used for the fabrication of field-effect transistors
with a nanowire channel [22, 23], which became the basis of biosensors for
detecting molecules and viruses [24], as well as for nanoscale field charge
sensors of an atomic force microscope for monitoring charge dynamics in
various structures [25].

In this Letter, we briefly describe a technology for fabricating
nanomechanical resonators from silicon-on-insulator based on the standard
microelectronics processes, only. The main parameters of the resonators, viz.,
its eigenfrequencies and quality factors were determined experimentally. The
dynamics of resonators in the linear and nonlinear regimes was investigated.

##  Oscillation properties of nanowires

The design of the most common types of nanomechanical resonators is based on a
suspended beam clamped at one or both ends. Such resonators have different
vibration modes of bending, twisting, deformation, etc. The bending mode is
the most interesting, since it produces the maximum response to an external
drive. This mode is easier to excite and detect by converting mechanical
vibrations into an electrical signal. The oscillation dynamics of a suspended
doubly clamped beam is generally described by the Euler-Bernoulli theory and
well approximated by the motion of a simple harmonic oscillator with small
damping [26]. For the fundamental mode, the equation of motion of the midpoint
of the beam excited by an external force has the following form:

$\ddot{x}+\gamma\dot{x}+\omega_{0}^{2}x=f(t),$ (2)

where $x$ is the displacement of the beam midpoint from the equilibrium
position (in the absence of an external force), $\gamma$ is the resonator
damping constant, $\omega_{0}=2\pi f_{0}$ is the angular resonance frequency,
and $f(t)$ is the external force.

For oscillations parallel to the plane of the substrate, the resonance
frequency of the fundamental bending mode is expressed by the following
expression [26]:

$f_{0}=1.03\sqrt{\frac{E}{\rho}}\frac{w}{L^{2}},$ (3)

where $E$ and $\rho$ are the Young’s modulus and density of the substrate
material, respectively, and $w$ and $L$ are the width and length of the
suspended nanowire. This equation is valid in the absence of the nanowire
tension, which in principle can arise when the sample is cooled down because
of the difference of the coefficients of thermal expansion of the materials
used.

The resonator remains in the linear regime at small mechanical displacements,
however, at a higher driving force, the resonator can enter the nonlinear
regime, which is accounted for by adding an extra term $\propto x^{3}$ to the
left-hand side of Eq. (2). In the nonlinear regime, the response curve becomes
asymmetric [27, 28, 29, 30]. Depending on the configuration of the resonant
system, the resonance frequency may increase (the resonator becomes “harder”)
or decrease (the resonator becomes “softer”), which is taken into account by
the sign of the cubic term. The resonator becomes a nonlinear system at
oscillation amplitudes exceeding the critical amplitude $a_{c}$ expressed as
[27]:

$a_{c}=2f_{0}\frac{L^{2}}{\pi}\sqrt{\frac{\rho\sqrt{3}}{EQ}},$ (4)

where $Q$ is the quality factor of the resonator.

To study such resonators, the magnetomotive method is commonly employed [20,
31]: the resonator is placed in a constant uniform magnetic field directed
perpendicular to the main axis of the nanowire (to the substrate) and an RF
current is passed through the nanowire. As a result, under the Lorentz force,
the nanowire starts to bend perpendicular to the direction of the magnetic
field and the direction of the current flow. The resonance frequencies of the
system can be found from the frequency dependence of the transmission
coefficient of the RF signal. Only odd modes can be detected by this method
since the induced emf vanishes for even modes due to the symmetry of the
system. The fundamental bending mode corresponding to the half of the
oscillation wavelength has the maximum displacement amplitude hence the
largest induced response.

## Sample fabrication

Fig. 1:  (color on-line) Fabrication stages of a silicon nanomechanical
resonator: (a) aluminum mask formed on top of the SOI wafer; (b) nanowire
structure in the top silicon layer after reactive ion etching; (c) suspended
structure of the nanowire after removal of the silicon oxide layer; (d) final
structure after metallization with a $20\,\mathrm{nm}$-thick layer of
aluminum.

Conventional SOI wafers with a thickness of the upper silicon layer of
$110\,\mathrm{nm}$, separated from the silicon substrate by a layer of the
$200\,\mathrm{nm}$-thick silicon oxide were used for the sample fabrication.
The fabrication process is similar to that used for sensors based on the
field-effect transistors with a nanowire channel [23, 32, 25] and involved
three stages of electron-beam lithography with a positive resist and also
reactive ion and wet etchings. The main fabrication stages are shown in Fig.
1.

Four resonant structures of different lengths connected in parallel were
formed on the chip. With this layout, it is possible to measure all resonators
in one cooldown, since the resonators do not affect each other due to the
significant difference in their resonance frequencies. The disadvantage of
this layout is that the incoming RF current is distributed among all
resonators, which does not allow to estimate of the absolute value of the
magnitude of the external force. A scanning electron microscope image of the
final structure with the designed lengths of $1$, $2$, $3$ and $5\,\mathrm{\mu
m}$, is shown in Fig. 2.

Fig. 2: (Color online) False color image of a structure of four nanowires in a
scanning electron microscope. Blue color denotes the silicon support
substrate; magenta – an insulating layer of silicon oxide; light-yellow –
silicon nanowires and wiring electrodes covered with a layer of aluminum, as
well as parasitic shadows on the substrate formed after aluminum deposition at
three different angles. The estimated resonators dimensions are (length
$\times$ width) $5\,\mathrm{\mu m}\times 115\,\mathrm{nm}$, $3\,\mathrm{\mu
m}\times 90\,\mathrm{nm}$, $2\,\mathrm{\mu m}\times 85\,\mathrm{nm}$,
$1\,\mathrm{\mu m}\times 80\,\mathrm{nm}$.

## Experimental results and discussion

Resonance characteristics of the fabricated structures were investigated in a
vacuum, in a magnetic field up to $5\,\mathrm{T}$ at a temperature of
$20\,\mathrm{mK}$. The used frequency range was up to $500\,\mathrm{MHz}$ with
a $200\,\mathrm{Hz}$ measurement bandwidth. RF signal from the network vector
analyzer was attenuated by $35\,\mathrm{dB}$ at cryogenic temperatures and
applied to the structure. The transmitted signal was amplified by
$40\,\mathrm{dB}$ at room temperature and detected. The power values in the
article correspond to the power applied to the measured structure. The shape
of the frequency characteristic of the signal transmission coefficient is
caused by the negative emf during the motion of a conductor in a magnetic
field, has a characteristic minimum at the resonant frequency as it shown in
Fig. 3. The symmetry of this dependence indicates the linear regime of the
resonator under the small external driving force. By approximating the
experimental data by the Lorentz function, the resonance frequency and
$Q$-factor of the resonators are determined. The measured values of the
resonance frequencies for the $5$, $3$ and $2\,\mathrm{\mu m}$-long nanowires
were found to be $32.46$, $71.99$ and $150.25\,\mathrm{MHz}$, respectively.
This agrees within the 2% error with the theoretical estimates obtained from
Eq. (3). The estimated value of the resonance frequency for the
$1\,\mathrm{\mu m}$-long nanowire is outside of the measured frequency range.

Fig. 3: (Color online) Dependence of the signal transmission coefficient of
the RF signal on the frequency for a suspended $5\,\mathrm{\mu m}$ nanowire
(black) and the approximation of experimental data by the Lorentz function
(red).

The mass of the resonator $\sim 1.2\times 10^{-13}\,\mathrm{g}$ is estimated
from its geometric dimensions. The minimum detectable frequency noise $\sim
8\,\mathrm{Hz}\,\mathrm{Hz}^{-1/2}$ is determined from the magnitude of the
signal transmission coefficient fluctuations
$10^{-3}\,\mathrm{a.\,u.}\,\mathrm{Hz}^{-1/2}$ and the greatest value of the
response gradient $125\times 10^{-6}\,\mathrm{a.\,u.}\,\mathrm{Hz^{-1}}$. By
plugging these numbers into Eq. (1) one can estimate the mass sensitivity in
the linear regime $\sim 6\times 10^{-20}\,\mathrm{g}\,\mathrm{Hz}^{-1/2}$,
which is comparable with the values of mass sensitivity obtained by other
groups [19]. The minimum detectable mass is $\sim 8.5\times
10^{-19}\,\mathrm{g}$, which corresponds, for example, to a silicon sphere of
radius $4.4\,\mathrm{nm}$. When the sensor operates in the nonlinear regime,
its mass sensitivity can be significantly improved [33].

Fig. 4: (a) Dependence of the oscillation amplitude of $5\,\mathrm{\mu m}$ on
the square of the magnetic field strength and its linear approximation
(straight line). (b) Dependence of the inverse $Q$-factor values of the
resonator on the square of the magnetic field strength and its linear
approximation.

As shown in [34], when a mechanical resonator is excited by a magnetomotive
method in the linear regime, the displacement amplitude is proportional to the
square of the magnetic field and to the amplitude of the RF current, $I\propto
P^{1/2}$, where $P$ is the RF power. Figure 4(a) shows the dependence of the
oscillation amplitude of the $5\,\mathrm{\mu m}$-long resonator on the square
of the magnetic field $B$. In magnetic fields up to $3.5\,\mathrm{T}$ the
experimental data are well approximated by the expected linear dependence of
the amplitude of mechanical oscillations on the magnitude of the magnetic
induction. With an increase in the magnetic field above $3.5\,\mathrm{T}$, the
resonator begins to enter the nonlinear regime and the amplitude increase
becomes weaker.

The width of the resonant response $\Delta f$ is determined by the losses both
in the resonator itself – internal losses, and in the measurement system –
external losses. For strong magnetic fields, the main ones are magnetomotive
losses which are proportional to the square of the magnetic field,
$Q^{-1}\propto\Delta f\propto B^{2}$. The measured quality factor increases
from $6.0\times 10^{3}$ at $5\,\mathrm{T}$ to $2.95\times 10^{4}$ at
$0.5\,\mathrm{T}$. The experimental dependence of $Q^{-1}$ on the square of
the magnitude of the magnetic field is shown in Fig. 4(b). The linearity of
this dependence indicates the dominance of magnetomotive losses over internal
ones. Such dependence is preserved throughout the investigated range of
magnetic fields. The internal $Q$-factor of the resonator is estimated to be
$3.62\times 10^{4}$ by the extrapolation to zero magnetic field. Such
$Q$-factor values exceed those obtained earlier in similar structures made
from SOI [35] and are comparable to the values in structures made of silicon
[36, 37], silicon coated with aluminum [38] and diamond [39]. It was confirmed
in the earlier experiments [36, 38, 39] that losses in resonators (which are
proportional to the inverse $Q$-factor values) decrease at temperatures below
$1\,\mathrm{K}$, but the nature of this dependence has not been understood
yet. The development of a reliable physical model describing losses in such
systems at low temperatures is at present the subject of scientific research.
In our experiments, measurements of the temperature dependence of the
$Q$-factor were not carried out, but it is obvious that the low temperature is
crucial for obtaining the high $Q$-factor of nanomechanical resonators.

Fig. 5: Dependence of the oscillation amplitude of the $2\,\mathrm{\mu
m}$-long resonator on the square root of the input power. The inset shows the
linear part of the entire dependence at low RF signal powers up to
$0.1\,\mathrm{a.\,u.}$. Straight line is the linear approximation of
experimental data.

Figure 5 shows the dependence of the amplitude of mechanical oscillations on
the square root of the power of the applied RF signal. The inset illustrates
the fact that at low powers the amplitude of the oscillations is proportional
to the oscillating force and the dynamics of the resonator is described by the
linear equation of the harmonic oscillator, Eq. (2). A slower increase in the
amplitude of the oscillations with increasing driving force
($P^{1/2}>0.1\,\mathrm{a.\,u.}$) indicates that the resonator becomes
“harder”. This also affects the shape of the response curve, which bends
towards the higher frequencies, as shown in Fig. 6. It is easy to show that
adding a cubic term $\propto x^{3}$ to Eq. (2) leads to a quadratic dependence
of the resonance frequency on the amplitude of the mechanical displacement
[27], which is observed in our experiments. The transition to the nonlinear
regime corresponds to the critical amplitude $a_{c}\sim 0.9\,\mathrm{nm}$
calculated from Eq. (4). From this, the largest amplitude of the resonator
oscillations $\sim 2.5\,\mathrm{nm}$ at the maximum pumping power of [28] is
found.

##  Conclusion

This Letter reports the resonance properties of the doubly clamped silicon
nanowires. By using the magnetomotive measurements scheme, we determined the
resonance frequencies of the fundamental bending modes of the nanowires in the
range of $30-150\,\mathrm{MHz}$. The intrinsic $Q$-factor of the
$5\,\mathrm{\mu m}$-long resonator was found to be $3.62\times 10^{4}$, which
is record high for resonators of similar dimensions fabricated from silicon-
on-insulator. The measured values of the resonance frequencies are in good
agreement with the estimates obtained from the Euler-Bernoulli theory. The
dynamics of the resonators is described qualitatively by the model of the
simple harmonic oscillator with a small dissipation and nonlinearity under the
external driving force. The fabrication method for resonant structures uses
standard silicon technological recipes only and can be used to fabricate high-
sensitivity mass, force and displacement sensors. The mass sensitivity of the
suspended $5\,\mathrm{\mu m}$-long nanowire is estimated to be $\sim 6\times
10^{-20}\,\mathrm{g}\,\mathrm{Hz}^{-1/2}$.

Fig. 6: (Color online) A family of response curves of the $2\,\mathrm{\mu
m}$-long resonator at high drive powers. When the RF power is increased, the
resonator becomes “harder”. The resonance frequency increases with the
amplitude of the resonator oscillations. The position of the maxima of the
resonance curves is approximated by a parabola, “backbone curve”, from which
the frequency of oscillations of the resonator in the linear regime is
determined.

The work was supported by the Russian Foundation for Basic Research, RFBR
Grant No. 16-29-03266 and the Royal Society (grant IEC\R3\170029). The
equipment of the “Educational Center of Lithography and Microscopy,” M. V.
Lomonosov Moscow State University was used.

## References

  * [1] Ekinci K. L., and Roukes M. L., Rev. Sci. Inst. 76, 061101, (2005).
  * [2] Zhang W-M., et al., Sensors 15, 26478, (2015).
  * [3] Arash B., Jiang J-W., and Rabczuk T., Appl. Phys. Rev. 2, 021301, (2015).
  * [4] Greenberg Ya. S., Pashkin Yu. A., and Il ichev E., Phys. Usp. 55 382, (2012).
  * [5] Shorokhov V. V., et al., Nanoscale 9, 613, (2017).
  * [6] Lovat G., et al., Nat. Nanotech. 12, 1050, (2017).
  * [7] Soldatov E. S., et al., Phys. Usp. 41 202, (1998)
  * [8] Bartsch S. T., Arp M., and Ionescu A. M., IEEE J. Elect. Dev. Soc. 2, 8, (2014).
  * [9] Ilic B., Yang Y., and Craighead H. G., Appl. Phys. Lett. 85, 2604, (2004).
  * [10] Mamin, H. and Rugar, D., Appl. Phys. Lett. 79, 3358, (2001).
  * [11] Zhao X., et al., Optics Express 20, 8535, (2012).
  * [12] Knobel R. G., and Cleland A. N., Nature 424, 291, (2003).
  * [13] Shevyrin A. A., et al., Appl. Phys. Lett. 106, 183110, (2015).
  * [14] Naik A., et al., Nature 443, 193, (2006).
  * [15] Teufel J., et al., Nature 475, 359, (2011).
  * [16] Harrabi K., et al., Appl. Phys. A 108, 7, (2012)
  * [17] Bradley D. I., et al., Sci. Rep. 7, 4876, (2017).
  * [18] Carr D. W., et al., Appl. Phys. Lett. 75, 920, (1999).
  * [19] Ekinci K. L., Huang X. M. H., and Roukes M. L, Appl. Phys. Lett. 84, 4469, (2004).
  * [20] Cleland A. N., and Roukes M. L., Appl. Phys. Lett. 69, 2653, (1996).
  * [21] Mori K., Silicon-On-Insulator (SOI) Technology, Chapter 14, 435, Woodhead Publishing, (2014).
  * [22] Presnov D. E., Amilonov S. V., and Krupenin V. A., Russian microelectronics, 41, 310, (2012)
  * [23] Presnov D. E., et al.,  Beilstein J. Nanotechnol. 4, 330, (2013).
  * [24] Rubtsova M., et al., “Biosensors 2016”, Procedia Technology, 27, 234, (2017).
  * [25] Trifonov A. S., et al., Ultramicroscopy 179, 33, (2017).
  * [26] Cleland A. N., Foundations of nanomechanics, Springer, Berlin, (2003)
  * [27] Nayfeh A. H. and Mook D. T.,  Nonlinear Oscillations., WILEY–VCH Verlag GmbH & Co. KGaA. (2007).
  * [28] Postma H. W. Ch., et al., Appl. Phys. Lett. 86, 223105, (2005).
  * [29] Tajaddodianfar F., Yazdi M. R. H., and Pishkenari H. N., Microsyst. Technol. 23, 1913, (2017).
  * [30] Laurent L., et al., Sens. and Act. A: Phys. 263, 326, (2017).
  * [31] Li T. F., et al., Appl. Phys. Lett. 92, 043112, (2008)
  * [32] Presnov D. E., Amitonov S. V. and Krupenin V. A., Russian Microelectronics 41, 364, (2012).
  * [33] Buks E., and Yurke B., Phys. Rev. E 74, 046619, (2017).
  * [34] Cleland A. N., and Roukes M. L., Sens. and Act., 72, 256, (1999).
  * [35] Yu L., et al., IEEE Trans. Nanotech. 11, 1093, (2012).
  * [36] Zolfagharkhani G., et al., Phys. Rev. B 72, 224101, (2005).
  * [37] Sulkko J., et al., Nano Lett. 10, 4884, (2010).
  * [38] Hoehne F., et al., Phys. Rev. B 81, 184112, (2010).
  * [39] Imboden M., Mohanty P., Phys. Rev. B 79, 125424, (2009).


Evolutionary status of W Vir pulsating variables

Yu. A. Fadeyev***E–mail<EMAIL_ADDRESS>
Institute of Astronomy, Russian Academy of Sciences, Pyatnitskaya ul. 48,
Moscow, 119017 Russia

Received October 12, 2020; revised October 16, 2020; accepted October 27, 2020

Abstract — Stellar evolution calculations for population II stars with initial
composition $Y_{0}=0.25$, $Z_{0}=10^{-3}$ and the initial stellar mass
$M_{0}=0.82M_{\odot}$ were carried out from the main sequence to the white
dwarf stage. Twelve AGB and post–AGB evolutionary sequences were computed with
different values of the parameter in the Blöcker mass loss rate formula
($0.01\leq\eta_{\mathrm{B}}\leq 0.12$). Selected models of evolutionary
sequences with masses $M=0.536M_{\odot}$, $0.530M_{\odot}$ and
$0.526M_{\odot}$ that experience the loop in the Hertzsprung–Russel diagram
due to the final helium flash were used as initial conditions for solution of
the equations of hydrodynamics describing radial stellar oscillations. The
region of instability to radial fundamental mode pulsations is shown to extend
from the asymptotic giant branch to effective temperatures as high as
$T_{\mathrm{eff}}\approx 6\times 10^{3}$ K. Pulsation periods of hydrodynamic
models are in the range from 15 to 50 day and agree with periods of W Vir
pulsating stars. The models of intermediate spectral type fundamental mode
pulsators with periods $\Pi>50$ day locate in the upper part of the
Hertzsprung–Russel diagram in the region of semiregular pulsating variables.
We conclude that W Vir pulsating variables are the low–mass post–AGB stars
that experience the final helium flash.

Keywords: stellar evolution; stellar pulsation; stars: variable and peculiar;
population II cepheids

## introduction

W Vir pulsating variables belong to population II cepheids and are observed
among field stars as well as in globular clusters. Periods of these variables
range from 8 to 35 day (Samus’ et al. 2017). W Vir stars are remarkable due to
strong hydrogen emission lines observed during rising light (Joy 1937).
Doubled metallic absorption lines detected simultaneously with hydrogen
emission are explained by the model of periodic shock wave propagating through
outer expanding layers of the stellar atmosphere (Abt 1954; Whitney 1956;
Wallerstein 1959). Periodic shock waves in the stellar atmosphere are due to
nonlinear effects in the large amplitude radial oscillations. Stronger
hydrogen emission is observed in RV Tau stars with higher luminosities and
longer periods. RV Tau semiregular pulsating variables are believed to be
connected with W Vir stars (Wallerstein 2002).

W Vir variability is thought to appear in low–metallicity low–mass double
shell source stars (Wallerstein and Cox 1984). Stellar evolution calculations
by Schwarzschild and Härm (1970) showed that thermal instability of the helium
burning shell may lead to occurence of the loop in the evolutionary track so
that the star crosses the cepheid instability strip in the Hertzsprung–Russel
(HR) diagram. More detailed evolutionary calculations by Gingold (1974, 1976)
and Sweigart et al. (1976) revealed that the necessary condition for occurence
of the blue loop in the HR diagram is that the helium flash occurs at the
onset of the post–AGB phase when the hydrogen envelope mass becomes as small
as $M_{\mathrm{env}}\approx 0.02M$, where $M$ is the mass of the star.

All studies of the evolutionary status of W Vir stars were done on the basis
of stellar evolution calculations whereas occurence of radial pulsations was
supposed from consideration of the loop on the evolutionary track with respect
to the instability strip of Cepheids and RR Lyr variables. For example,
Gingold (1976) extrapolated results of RR Lyr linear model calculations by
Tuggle and Iben (1972) to higher luminosities. On the other hand, computations
of nonlinear pulsation models of W Vir stars (Fadeyev and Fokin 1985; Buchler
and Kovacs 1987, Kovacs and Buchler 1988; Fadeyev and Muthsam 1990; Serre et
al. 1996; Smolec 2016) were carried out without a direct relation between the
stellar evolution models and the initial conditions of radiation hydrodynamics
calculations. Below we show that pulsations of W Vir stars strongly depend on
the mass and luminosity because of the low gas density in the envelope of the
pulsating star and significant contribution of radiation energy density to the
internal energy of stellar matter. For example, variation of the stellar mass
within 4% is accompanied by change of the luminosity by a factor of two.

The goal of the present study is to establish evolutionary status of W Vir
pulsating variables on the basis of consistent stellar evolution and nonlinear
stellar pulsation calculations. Methods of the solution are discussed in our
earlier papers devoted to studies of classical Cepheids (Fadeyev 2018) and
post–AGB stars (Fadeyev 2019a, 2019b).

## results of stellar evolution calculations

In the present work we consider stellar evolution from the main sequence to
the stage of a white dwarf with luminosity $L\sim L_{\odot}$. Initial
mass–fraction abundances of helium and heavier elements were assumed to agree
with those of typical globular clusters: $Y_{0}=0.25$ (Salaris et al. 2004)
and $Z_{0}=0.001$ (Catelan 2004). The initial mass of the star on the zero age
main sequence is $M_{0}=0.82M_{\odot}$ so that the age of W Vir models is
$t_{\mathrm{ev}}=1.244\times 10^{10}$ yr.

Evolutionary sequences were computed with the program MESA version 12778
(Paxton et al. 2019). In order to calculate the thermonuclear energy release
and evolution of isotope abyndances we employed the reaction network
’pp_cno_extras_o18_ne22.net’ with 26 nuclei from hydrogen ${}^{1}\mathrm{H}$
to magnesium ${}^{24}\mathrm{Mg}$ coupled by 81 reactions. The reaction rates
were calculated with the JINA Reaclib data (Cyburt et al. 2010). Convective
mixing was treated according to the theory of Böhm–Vitense (1958) for the
mixing length to pressure scale height ratio
$\alpha_{\mathrm{MLT}}=\Lambda/H_{\mathrm{P}}=2.0$. Mixing due to overshooting
at the boundaries of convective stability was computed according to Herwig
(2000) with parameters $f=0.016$ and $f_{0}=0.004$. Evolutionary calculations
of the core helium burning stage with option ’conv_premix_avoid_increase’
allowed us to significantly diminish the role of breathing pulses and to avoid
appearence of spurious loops on the evolutionary track. The mass loss rate at
evolutionary stages preceding AGB was evaluated by the Reimers (1975) formula
with parameter $\eta_{\mathrm{R}}=0.5$, whereas at later stages we used the
mass loss rate formula by Blöcker (1995). Bearing in mind existing
uncertainties in mass loss rates of AGB stars we computed twelve AGB and
post–AGB evolutionary sequences with parameter $\eta_{\mathrm{B}}$
equidistantly distributed in the range $0.01\leq\eta_{\mathrm{B}}\leq 0.12$.

Results of stellar evolution calculations are illustrated in Fig. 1, where
three tracks of AGB and post–AGB stages are shown in the HR diagram for the
mass loss parameters $\eta_{\mathrm{B}}=0.03$, 0.07 and 0.1. According to
Miller Bertolami (2016) we assume that the post–AGB stage begins when the mass
of the hydrogen envelope is $M_{\mathrm{env}}\leq 0.01M$. Transition from the
AGB to post–AGB stage is marked on each track by the open circle, whereas the
filled circle corresponds to the peak of the helium–flash luminosity.

Of most interest for the present study are the evolutionary changes of the
star in the effective temperature range $4\times 10^{3}\>\mathrm{K}\lesssim
T_{\mathrm{eff}}\lesssim 6\times 10^{3}\>\mathrm{K}$. Evolutionary tracks
within this effective temperature interval are shown with enlarged scale in
Fig. 2.

Main properties of evolutionary sequences are listed in Table 1, where
$n_{\mathrm{TP}}$ and $\langle\Delta t_{ip}\rangle$ are the number of thermal
pulses during the AGB stage and the mean time interval between two thermal
pulses, $M$ and $L$ are the mass and luminosity of the star at the onset of
the post–AGB stage when the mass of the hydrogen envelope is
$M_{\mathrm{env}}=0.01M$, $M_{\mathrm{env}}^{\star}$ is the mass of the
hydrogen envelope at the peak of the helium luminosity during the final helium
flash. The sign ’$-$’ in the last column implies absence of the blue loop due
to the final helium flash in the HR diagram.

Increase of the mass loss rate leads to the shorter AGB lifetime and to the
smaller number of thermal pulses $n_{\mathrm{TP}}$. Two temporal dependences
of the helium–burning shell luminosity $L_{\mathrm{He}}$ for evolutionary
sequences $\eta_{\mathrm{B}}=0.01$ и $\eta_{\mathrm{B}}=0.1$ are shown in Fig.
3, where for the sake of graphical representation the evolutionary time
$t_{\mathrm{ev}}$ is set to zero at the helium peak luminosity of the first
thermal pulse $t_{\mathrm{TP}1}$. As seen from these plots, the TP–AGB
lifetime reduces by nearly two times when the mass loss rate increases by ten
times although the mean time interval between flashes remains almost
unchanged.

The necessary condition for occurence of the loop in the evolutionary track
during the post–AGB stage is that the final thermal pulse occurs at the mass
of the hydrogen envelope $M_{\mathrm{env}}\approx 0.02M$ (Gingold 1974). As
seen in Table 1, results of our computations agree with this condition because
in the three evolutionary sequences without the loop on the evolutionary track
($\eta_{\mathrm{B}}=0.05$, 0.06 and 0.12) the mass of the hydrogen envelope
during the final thermal pulse is significantly larger. It should be noted
that for $\eta_{\mathrm{B}}=0.08$ the mass of the hydrogen envelope at the
final thermal pulse is $M_{\mathrm{env}}=0.020M$ but the loop of the
evolutionary track extends in the HR diagram to effective temperatures as high
as $T_{\mathrm{eff}}\approx 5000$ K. That is why in Table 1 this evolutionary
sequence is marked as that without the loop. This example allows us to
conclude that the sufficiently small value of the mass of the hydrogen
envelope during the final helium flash is not sufficient for occurence of the
extended loop on the evolutionary track.

## results of stellar pulsation calculations

In consistent stellar evolution and nonlinear stellar pulsation calculations
selected evolutionary models are used as initial conditions for solution of
the equations of radiation hydrodynamics. In the present work the system of
the equations of hydrodynamics is extended due to transport equations of
time–dependent turbulent convection (Kuhfuß 1986) and is discussed in our
earlier paper (Fadeyev 2013).

Solution of the Cauchy problem for the equations of hydrodynamics implies that
the role of initial perturbations is played by interpolation errors that
appear when the physical variables of the evolutionary model (the radius, the
luminosity, the temperature and the pressure) are used for calculation in the
Lagrangian grid of the hydrodynamic model. Evolutionary calculations of the
final AGB stage were done with the number of mass zones $N_{\mathrm{MESA}}\sim
10^{4}$, whereas the number of mass zones of the hydrodynamic models was
$N=500$. The inner boundary of the hydrodynamic model was determined in the
layers with radiative transfer and the gas temperature $T\sim 10^{6}$ K. The
radius of the inner boundary is $r\lesssim 0.05R$, where $R$ is the radius of
the evolutionary model. The radius and the luminosity of the inner boundary
were assumed to be constant on the time interval of the Cauchy problem
solution.

The solution of the equations of hydrodynamics describes decaying oscillations
if the star model is stable to radial pulsations but in the case of the
pulsational instability the solution describes oscillations with exponentially
increasing amplitude until the limiting amplitude is reached. The main
difference between classical Cepheids and W Vir pulsating variables is that
the stars of the latter group have significantly higher rates of the amplitude
growth. For example, typical growth rate of the kinetic energy of pulsation
motions in W Vir models is $\eta=\Pi d\ln E_{\mathrm{K,max}}/dt\approx 1$,
where $\Pi$ is the pulsation period and $E_{\mathrm{K,max}}$ is the maximum
kinetic energy reached during the period.

Hydrodynamic calculations of nonlinear stellar pulsations were carried out for
selected models of evolutionary sequences computed with mass loss rate
parameters $\eta_{\mathrm{B}}=0.03$, 0.07 and 0.1. Instability to radial
fundamental mode pulsations was found for stars between AGB and effective
temperatures as high as $T_{\mathrm{eff}}\approx 6000$ K. Estimates of
$T_{\mathrm{eff}}$ corresponding to the blue edge of pulsational instability
($\eta=0$) are shown by filled squares in Figs. 1 and 2.

In order to establish the evolutionary status of W Vir stars we have to
compare the pulsation periods of hydrodynamic models with periods of observed
light variations. As seen in Fig. 1, for the parameter
$\eta_{\mathrm{B}}=0.03$ the evolutionary track crosses the instability edge
only once during the evolution from AGB. The extension of the red loop in the
HR diagram is not sufficient and the track turns back at the effective
temperature $T_{\mathrm{eff}}\approx 7300$ K. Hydrodynamic models of
evolutionary sequences $\eta_{\mathrm{B}}=0.07$ and $\eta_{\mathrm{B}}=0.1$
appear in the region of pulsational instability two times: during evolution
from AGB and when the evolutionary track loops to the red giant region in the
HR diagram and the effective temperature becomes $T_{\mathrm{eff}}<6000$ K.

The plots of the period change rate as a function of period for the stage of
evolution from AGB to the first crossing of the blue edge of the instability
region are shown in Fig. 4. The initial point of each plot represents the
hydrodynamic model with effective temperature $T_{\mathrm{eff}}\approx 5000$ K
($\eta_{\mathrm{B}}=0.03$, 0.07) and $T_{\mathrm{eff}}\approx 4500$ K
($\eta_{\mathrm{B}}=0.1$). Evolution proceeds with decreasing period from left
to right to the blue edge of the pulsation instability region and lasts from
2200 yr for $\eta_{\mathrm{B}}=0.03$ to 3200 yr for $\eta_{\mathrm{B}}=0.1$.

Variations of $\Pi$ and $\dot{\Pi}$ between the second crossing of the
instability edge and the maximum period are shown in Fig. 5a. In models of the
evolutionary sequence $\eta_{\mathrm{B}}=0.07$ the pulsation period increases
from $\Pi=17.6$ day to $\Pi=21.2$ day during $\approx 5800$ yr and the mass of
the hydrogen envelope ($M_{\mathrm{env}}=0.011M$) almost does not change. In
hydrodynamic models of the evolutionary sequence $\eta_{\mathrm{B}}=0.1$ the
period of radial oscillations increases from 28.3 day to 89 day during 165 yr,
whereas the mass of the hydrogen envelope decreases from
$M_{\mathrm{env}}=4.3\times 10^{-3}$ to $M_{\mathrm{env}}\approx 10^{-3}$.
Thus the W Vir star with mass $M=0.526M_{\odot}$ becomes the semiregular
pulsating variable during one and a half century.

Hydrodynamic calculations of radial pulsations for the stage with decreasing
period after the second crossing of the blue edge of the pulsation instability
region were done only for the models of the evolutionary sequence
$\eta_{\mathrm{B}}=0.07$ because for $\eta_{\mathrm{B}}=0.1$ pulsation periods
significantly exceed those of W Vir pulsating variables. Results of these
calculations are shown in Fig. 5b. The mass of the hydrogen envelope reduces
from $M_{\mathrm{env}}=0.011M$ to $4.6\times 10^{-3}M$ whereas the pulsation
period decreases from 21.4 day to 12.7 day during $\approx 9\times 10^{3}$ yr.

Nonlinearity of radial pulsations in W Vir variables is responsible for the
absence of strict repetition of light variations which is typical for the most
of classical Cepheids. This feature is a substantial obstacle in determination
of the rates of period change so that only for a few W Vir stars the
successfull measurements of $\dot{\Pi}$ have been done by now. To compare
results of our calculations with observations we considered two W Vir
pulsating variables with most reliable $O-C$ diagrams. One of them is CC Lyr
with period $\Pi=24.01$ day and the rate of period increase $\dot{\Pi}=229.3$
s/yr (Berdnikov et al. 2020) and another is V1303 Cyg with period $\Pi=18.45$
day and the rate of period decrease $\dot{\Pi}=-40$ s/yr (Barton 1986). As
seen in Fig. 5, results of our calculations do not contradict the
observations.

## conclusions

Presented above results of consistent stellar evolution and nonlinear stellar
pulsation calculations allow us conclude that W Vir pulsating variables are
the post–AGB stars that experience the final helium flash. Only in this case
the luminosity of the post–AGB star is in the range $10^{3}L_{\odot}\lesssim
L\lesssim 2\times 10^{3}L_{\odot}$, whereas the periods of radial oscillations
($10\leavevmode\nobreak\ \textrm{day}\lesssim\Pi\lesssim
50\leavevmode\nobreak\ \textrm{day}$) correspond to periods of light
variations of W Vir stars. It should be noted that results of our calculations
confirm the conclusion of Gingold (1974) that the helium shell flash before
the onset of the post–AGB stage is the necessary condition for occurence of W
Vir variability. The term ’final thermal pulse’ appeared later in papers by
Iben (1982) and Schoenberner (1983).

Evolution of the post–AGB star after the final helium flash proceeds in the
thermal timescale of the hydrogen envelope. The small envelope mass and large
contribution of the radiation field in the internal energy of the stellar
matter are responsible for rapid evolution. Most demonstrative illustration of
rapid evolutionary changes is given by hydrodynamic models of the evolutionary
sequence $\eta_{\mathrm{B}}=0.1$ after the second crossing of the pulsation
instability edge. One of such objects with period change in time scale of
$\sim 10^{2}$ yr is V725 Sgr. In the beginning of the twentieth century this
pulsating variable was known as the population II cepheid but later this star
became the semiregular pulsating variable with period $80\leavevmode\nobreak\
\textrm{day}<\Pi<100\leavevmode\nobreak\ \textrm{day}$ (Percy 2006).

## references

  1. 1.

H.A. Abt, Astrophys. J. Suppl. Ser. 1, 63 (1954).

  2. 2.

A.S. Barton, J. Am. Associat. Var. Star Observ. 15, 246 (1986).

  3. 3.

L.N. Berdnikov, A.M. Yakob, E.N. Pastukhova, Astron. Lett. 46, 630 (2020).

  4. 4.

T. Blöcker, Astron. Astrophys. 297, 727 (1995).

  5. 5.

E. Böhm–Vitense, Zeitschrift für Astrophys. 46, 108 (1958).

  6. 6.

J.R. Buchler and G. Kovacs, Astrophys. J. 320, L57 (1987).

  7. 7.

M. Catelan, Astrophys. J. 600, 409 (2004).

  8. 8.

R.H. Cyburt, A.M. Amthor, R. Ferguson, Z. Meisel, K. Smith, S. Warren, A.
Heger, R.D. Hoffman, T. Rauscher, A. Sakharuk, H. Schatz, F.K. Thielemann, and
M. Wiescher, Astrophys. J. Suppl. Ser. 189, 240 (2010).

  9. 9.

Yu.A. Fadeyev, Astron. Lett. 39, 306 (2013).

  10. 10.

Yu.A. Fadeyev, Astron. Lett. 44, 782 (2018).

  11. 11.

Yu.A. Fadeyev, Astron. Lett. 45, 521 (2019a).

  12. 12.

Yu.A. Fadeyev, Astron. Lett. 45, 655 (2019b).

  13. 13.

Yu.A. Fadeyev and A.B. Fokin, Astrophys. Space Sci. 111, 355 (1985).

  14. 14.

Yu.A. Fadeyev and H. Muthsam, Astron. Astrophys. 234, 188 (1990).

  15. 15.

R.A. Gingold, Astrophys. J. 193, 177 (1974).

  16. 16.

R.A. Gingold, Astrophys. J. 204, 116 (1976).

  17. 17.

F. Herwig, Astron. Astrophys. 360, 952 (2000).

  18. 18.

I. Iben, 260, 821 (1982).

  19. 19.

A.H. Joy, Astrophys. J. 86, 363 (1937).

  20. 20.

G. Kovacs and J.R. Buchler, Astrophys. J. 334, 971 (1988).

  21. 21.

R. Kuhfuß, Astron. Astrophys. 160, 116 (1986).

  22. 22.

M.M. Miller Bertolami, Astron. Astrophys. 588, A25 (2016).

  23. 23.

B. Paxton, R. Smolec, J. Schwab, A. Gautschy, L. Bildsten, M. Cantiello, A.
Dotter, R. Farmer, J.A. Goldberg, A.S. Jermyn, S.M. Kanbur, P. Marchant, A.
Thoul, R.H.D. Townsend, W.M. Wolf, M. Zhang, and F.X. Timmes, Astrophys. J.
Suppl. Ser. 243, 10 (2019).

  24. 24.

J.R. Percy, A. Molak, H. Lund, D. Overbeek, A.F. Wehlau, and P.F. Williams,
Publ. Astron, Soc. Pacific 118, 805 (2006).

  25. 25.

D. Reimers, Problems in stellar atmospheres and envelopes (Ed. B. Baschek,
W.H. Kegel, G. Traving, New York: Springer-Verlag, 1975), p. 229.

  26. 26.

M. Salaris, M. Riello, S. Cassisi, and G. Piotto, Astron. Astrophys. 420, 911
(2004).

  27. 27.

N.N. Samus’, E.V. Kazarovets, O.V. Durlevich, N.N. Kireeva, and E.N.
Pastukhova, Astron. Rep. 61, 80 (2017).

  28. 28.

D. Schoenberner, 272, 708 (1983).

  29. 29.

M. Schwarzschild and R. Härm, Astrophys. J. 160, 341 (1970).

  30. 30.

T. Serre, Z. Kollath and J.R. Buchler, Astron. Astrophys. 311, 845 (1996).

  31. 31.

R. Smolec, MNRAS 456, 3475 (2016).

  32. 32.

A.V. Sweigart, J.G. Mengel, and P. Demarque, Astron. Astrophys. 30, 13 (1974).

  33. 33.

R.S. Tuggle and I. Iben, Astrophys. J. 178, 455 (1972).

  34. 34.

G. Wallerstein, Astrophys. J. 130, 560 (1959).

  35. 35.

G. Wallerstein, Publ. Astron, Soc. Pacific 114, 689 (2002).

  36. 36.

G. Wallerstein and A.N. Cox, Publ. Astron, Soc. Pacific 96, 677 (1984).

  37. 37.

C. Whitney, Annales d’Astrophysique 19, 142 (1956).

Table 1: Main properties of evolutionary sequences $\eta_{\mathrm{B}}$ | $n_{\mathrm{TP}}$ |  $\langle\Delta t_{\mathrm{ip}}\rangle,\ 10^{6}$ yr | $M/M_{\odot}$ | $L/L_{\odot}$ | $M_{\mathrm{env}}^{\star}/M_{\odot}$ |
---|---|---|---|---|---|---
0.01 | 7 | 0.234 | 0.550 | 3342 | 0.009 |
0.02 | 6 | 0.236 | 0.541 | 3093 | 0.005 |
0.03 | 5 | 0.237 | 0.536 | 1883 | 0.025 |
0.04 | 5 | 0.238 | 0.535 | 2818 | 0.004 |
0.05 | 4 | 0.235 | 0.533 | 2637 | 0.062 | $-$
0.06 | 3 | 0.376 | 0.531 | 2291 | 0.033 | $-$
0.07 | 4 | 0.233 | 0.529 | 1866 | 0.035 |
0.08 | 4 | 0.232 | 0.529 | 5448 | 0.020 | $-$
0.09 | 4 | 0.237 | 0.528 | 2372 | 0.005 |
0.10 | 4 | 0.236 | 0.528 | 2315 | 0.004 |
0.11 | 4 | 0.235 | 0.527 | 2218 | 0.004 |
0.12 | 3 | 0.233 | 0.526 | 2040 | 0.061 | $-$

Fig. 1: Evolutionary tracks in the HR diagram for the mass loss parameters
$\eta_{\mathrm{B}}=0.1$ (solid line), $\eta_{\mathrm{B}}=0.07$ (dashed line)
and $\eta_{\mathrm{B}}=0.03$ (dotted line). Filled circles mark maxima of the
helium burning shell during the final helium flash, open circles mark the
onset of the post–AGB stage, filled squares mark the blue edge of the
pulsation instability region determined from hydrodynamic computations.

Fig. 2: Same as Fig. 1 but with enlarged scales of the HR diagram area
confined by the asymptotic giant branch and the blue edge of the pulsation
instability region.

Fig. 3: The helium–burning shell luminosity $L_{\mathrm{He}}$ during the AGB
stage for the mass loss rate parameters $\eta_{\mathrm{B}}=0.01$ (a) and
$\eta_{\mathrm{B}}=0.1$ (b). Open circles on the plots mark the onset of the
post–AGB stage when $M_{\mathrm{env}}=0.01M$. The evolutionary time
$t_{\mathrm{ev}}$ on the horizontal axis is set to zero at the helium peak
luminosity of the first thermal pulse $t_{\mathrm{TP}1}$.

Fig. 4: The rate of period change $\dot{\Pi}$ as a function of period $\Pi$
for the evolutionary stage from AGB up to the first crossing of the blue edge
of the pulsation instability region.

Fig. 5: The rate of period change $\dot{\Pi}$ as a function of period $\Pi$
after the second crossing of the blue edge of the pulsation instability region
during the stage of period increase (a) and period decrease (b). Evolutionary
sequences $\eta_{\mathrm{B}}=0.07$ and $\eta_{\mathrm{B}}=0.1$ are shown by
solid and dashed lines. Observational estimates of $\Pi$ and $\dot{\Pi}$ for
CC Lyr (Berdnikov et al. 2020) and V1303 Sgr (Barton 1986) are shown by filled
circles.