# CEHR-GPT: Generating Electronic Health Records with Chronological Patient
Timelines

Chao Pang111indicates equal contribution Department of Biomedical Informatics,
Columbia University Irving Medical Center  Observational Health Data Sciences
and Informatics  Xinzhuo Jiang∗ Department of Biomedical Informatics, Columbia
University Irving Medical Center  Observational Health Data Sciences and
Informatics  Nishanth Parameshwar Pavinkurve  Department of Biomedical
Informatics, Columbia University Irving Medical Center  Observational Health
Data Sciences and Informatics  Krishna S. Kalluri  Department of Biomedical
Informatics, Columbia University Irving Medical Center  Observational Health
Data Sciences and Informatics  Elise L. Minto  Department of Biomedical
Informatics, Columbia University Irving Medical Center  Observational Health
Data Sciences and Informatics  Jason Patterson  Department of Biomedical
Informatics, Columbia University Irving Medical Center  Observational Health
Data Sciences and Informatics  Linying Zhang  Institute for Informatics, Data
Science, and Biostatistics. Washington University in St Louis  Observational
Health Data Sciences and Informatics  George Hripcsak  Department of
Biomedical Informatics, Columbia University Irving Medical Center  Medical
Informatics Services, New York-Presbyterian Hospital  Observational Health
Data Sciences and Informatics  Noémie Elhadad  Department of Biomedical
Informatics, Columbia University Irving Medical Center  Medical Informatics
Services, New York-Presbyterian Hospital  Observational Health Data Sciences
and Informatics  Karthik Natarajan  Department of Biomedical Informatics,
Columbia University Irving Medical Center  Medical Informatics Services, New
York-Presbyterian Hospital  Observational Health Data Sciences and Informatics

## Abstract

Synthetic Electronic Health Records (EHR) have emerged as a pivotal tool in
advancing healthcare applications and machine learning models, particularly
for researchers without direct access to healthcare data. Although existing
methods, like rule-based approaches and generative adversarial networks
(GANs), generate synthetic data that resembles real-world EHR data, these
methods often use a tabular format, disregarding temporal dependencies in
patient histories and limiting data replication. Recently, there has been a
growing interest in leveraging Generative Pre-trained Transformers (GPT) for
EHR data. This enables applications like disease progression analysis,
population estimation, counterfactual reasoning, and synthetic data
generation. In this work, we focus on synthetic data generation and
demonstrate the capability of training a GPT model using a particular patient
representation derived from CEHR-BERT, enabling us to generate patient
sequences that can be seamlessly converted to the Observational Medical
Outcomes Partnership (OMOP) data format.

## Keywords

Generative Pre-trained Transformer, Synthetic Electronic Health Records,
Patient Representation, Observational Medical Outcomes Partnership - Common
Data Model, Observational Health Data Sciences and Informatics

##  1 Introduction

Access to electronic health records (EHRs) is fundamental to healthcare
research, drug surveillance, clinical machine learning, and system
development. However, the use of real-world EHR data comes with considerable
challenges such as privacy and security issues, institutional consent, and
restrictions on data sharing. Synthetic data emerges as a promising solution,
offering a more expedient and secure pathway to healthcare information, which
could accelerate progress across various sectors, including academic research,
clinical settings, and the pharmaceutical industry [1].

Synthetic data is not real data, that is, it doesn’t relate to any specific
individual. However, it mimics the statistical characteristics and journeys of
specific patient populations. Synthetic data enables a broader range of
researchers to answer their questions of interest without going through the
cumbersome process of accessing real data and worrying about patient privacy
[2]. In recent years, many machine learning, specifically deep learning and
generative artificial intelligence (AI) models, have been developed to derive
synthetic data from real EHR data. [3] However, most existing methods for
synthetic EHR data generation fail to adequately capture the temporal
dependencies that are often critical in medical scenarios. These temporal
aspects, such as medication schedules, symptom progression, and lab result
timelines, are vital for a comprehensive understanding of patient health
trajectories and for developing effective treatment strategies. An ideal
synthetic dataset derived from institutional data should maintain the inherent
correlations among time-series features, thus enabling researchers to
externally validate machine learning models in different populations.
Crucially, the synthetic dataset must preserve accurate patient timelines, as
predictive tasks are highly susceptible to temporal variations. A synthetic
dataset is considered to exhibit comparable machine learning utility to the
original data if it meets two key criteria: 1) it demonstrates similar outcome
prevalence to the source data; 2) machine learning models trained on the
synthetic data achieve performance metrics akin to those trained with the
original data.

The majority of the existing research focuses on developing new deep learning
models in generative EHR research but without adequate emphasis on retaining
accurate temporal information [4, 5, 6]. Unfortunately, synthetic EHR datasets
developed as such will not support use cases that require the accurate
construction of a patient timeline e.g., 30-day readmission, one-year risk of
heart failure, and disease progression. This limits existing work to only
perform simple code prediction in their evaluations instead of comprehensive
phenotype predictions. Another challenge in using synthetic EHR data in
practice is its difficulty with dissemination due to a lack of standards.
Synthetic patient sequences cannot be widely adopted for analyses without use
of a common data model, however, none of the existing works have included such
a component in their frameworks to present the synthetic data in an easy-to-
consume fashion.

In our view, time-series synthetic data should not only capture the underlying
characteristics of heterogeneous EHRs but also satisfy the following temporal
requirements, 1) a matching distribution of the starting age; 2) a matching
distribution of the starting year; 3) a matching distribution of the inpatient
duration; 4) a matching distribution of time intervals between neighboring
visits. Furthermore, synthetic EHR data should be stored in common data models
such as the Observational Medical Outcomes Partnership (OMOP) Common Data
Model, which is used in many large data networks, [7] for easy dissemination
and consumption. Although creating such a time-series EHR dataset seems to be
a challenging task, we think that this problem can be solved through a patient
representation approach. The key is to focus on designing a good patient
representation rather than creating a sophisticated architecture to model time
and medical events simultaneously. In this paper, we present the CEHR-GPT
framework for building an end-to-end workflow to generate time-series
synthetic EHR data. Our contributions are summarized below,

  * •

We design a novel patient representation that not only keeps track of the
visit types and discharge facility for inpatient visits but also preserves all
temporal information including patients’ starting year, starting age, time
intervals between visits, and inpatient visit spans. To the best of our
knowledge, this is the first time that temporal information is fully retained.

  * •

We treat patient sequence generation as a language modeling problem, which
allowed us to use the state-of-the-art language model Generative Pre-trained
Transformers (GPT) to learn the distribution of patient sequences to generate
new synthetic sequences [8, 9].

  * •

We converted synthetic sequences to the common data format OMOP with almost no
loss of temporal information. Synthetic OMOP can be easily evaluated using the
OHDSI tools and disseminated to others.

  * •

We evaluated the synthetic EHR data on three levels, dimension-wise
distribution (marginal distribution), co-occurrence relationship, and machine
learning model performance metrics.

##  2 Related work

With the adoption of Generative Adversarial Networks (GANs) [10], researchers
have found creative ways to generate synthetic EHR data. Since 2017, several
groups have applied GANs to tabular EHRs and developed several evaluation and
privacy metrics to quantify the performance of GANs [3]. Despite the success,
one limitation is that the tabular format fails to capture the temporal nature
of EHR data because they are constructed from patient histories using a bag-
of-words approach. It was not until 2020 that researchers started developing
new GAN architectures to tackle the time series data. Dual adversarial
autoencoder (DAAE) [6] used a combination of a variational autoencoder (VAE)
and two GAN components, where the inner GAN was trained to replicate the
encoded representation generated by the encoder, and the outer GAN was trained
to generate realistic-looking patient sequences in addition to the
reconstruction error. Another model called EHR-M-GAN [5] employed a similar
autoencoder architecture with two main differences 1) they used a dual-VAE
framework to handle the continuous and discrete valued features, where the
continuous and discrete representations were generated by encoders first and
then collectively used for decoding; 2) in the GAN generator, they used two
parallel recurrent neural networks (one for sampling continuous noise and one
for discrete noise) with a so-called Bilateral LSTM cell to allow the
continuous and discrete noise vectors to interact with each other to generate
better sampling vectors. Although these GANs took the temporal order of events
into consideration, they did not generate timestamps for visits or medical
events, therefore limiting the use of the synthetic data.

Improving upon the previous works, a two-stage learning algorithm (dependency
learning and conditional simulation) named SynTEG was proposed to generate
timestamped synthetic data [11]. In dependency learning, transformer encoders
learned visit representations, which were used to feed into an recurrent
neural network (RNN) model to learn the dependencies between visits. Two self-
supervising learning tasks were used for training – prediction of the
timestamp and diagnoses at the next visit. In the conditional simulation, the
learned visit representations generated from the previous step were used to
train a conditional GAN to generate the diagnosis codes at each visit. This
approach achieved superior performance in the time-sensitive evaluations over
the previous methods. However, there are still ongoing challenges that have
not been addressed in their work, 1) other EHR data was not utilized because
only diagnosis codes were included in training; 2) the visits were assumed to
start and end on the same day, which would fail to model in-patient visits
that normally span days; therefore, the constructed patient timeline would be
inaccurate; 3) the synthetic data didn’t include visit types and discharge
facilities.

Until now, almost all the existing approaches used some variation of GAN for
learning the data distribution, but unfortunately, GANs are notoriously
difficult to train and easily subject to mode collapse. Despite the recent
advancements in optimization techniques such as Wasserstein-GAN [12], it would
require a significant amount of time to tune hyperparameters to properly train
GANs. As the previous works have demonstrated [13, 14, 15], patient sequence
generation could be conceptually represented as a language modeling problem.
Foresight [16] adapted GPT to forecast patient trajectories. Their methodused
a name entity recognition (NER) tool to extract medical concepts from
discharge summaries, based on which a patient sequence was constructed
chronologically, and then they trained a standard GPT model using all patient
sequences constructed from the previous step. To forecast future events, they
fed a patient history to prompt the model and employed a Monte Carlo sampling
strategy to calculate the probability of developing certain conditions. One
limitation of this work was that the model could not predict when a certain
condition would happen due to a lack of temporal information in their patient
sequence. Nevertheless, their work demonstrated the potential of using GPT for
modeling patient sequences. We aim to address this limitations in our work.

##  3 Methods

In Figure 1, we present a framework for generating synthetic EHR data from an
OMOP source. To retain the temporal dependencies, we opted to work directly
with time-series patient sequences instead of using a bag-of-words (BOW)
representation by building on our previous work [13]. We first encoded the
OMOP data into patient sequences using a specific patient representation,
described later in this section. Secondly, we trained a generative model on
the converted patient sequences and utilized it to generate new synthetic
patient sequences.

Figure 1: Overall architecture. The OMOP data is first converted to patient
sequences by an OMOP encoder based on the patient representation that
preserves demographics, visit types, and temporal intervals between visits.
Then a generative model is trained to learn the sequence distribution in order
to generate new sequences. Next, the generated sequences are converted back to
the OMOP format using an OMOP decoder.

Finally, we fed the synthetic patient sequences into an OMOP decoder to create
a synthetic OMOP dataset. Furthermore, an evaluation procedure was developed
to assess the similarity between the synthetic OMOP and the source OMOP data.

###  3.1 Patient Representation

We designed a patient representation in CEHR-BERT [13] that captures medically
relevant events and their timelines while exhibiting certain characteristics
of a sentence. In order to fully leverage Large Language Models (LLM) on
patient sequences, we further extended this patient representation to include
demographic information, patient history, and temporal dependencies as shown
in Figure 2. The start of the sequence defines the demographic prompt
containing EHR start year, age, gender, and race. It is followed by visit
blocks separated by artificial time tokens (ATT) representing time intervals
in days, e.g., $D_{1}$ represents an interval of 1 day. For time intervals
surpassing 1080 days, we grouped these into a single Long Term (LT) token, a
decision guided by the low occurrence rate in this time frame. Each visit
block starts with a visit type token (VTT) to signify the type of visit and is
then followed by domain records arranged in chronological order. In the case
of inpatient visits, a distinct inpatient ATT (IATT) was inserted between
neighboring inpatient spans, defined as the groups of records that occurred on
the same day. These are distinct from ATT and are used for capturing the time
between multiple events characterized by concepts within the same visit. In
addition, a discharge facility code (e.g. discharge home and long-term care)
was inserted at the end of the inpatient visit. The use of IATT, which is
distinct from ATT, was necessary since they are attributed to two different
contexts, and resulted in better performance than using the same ATT across
both contexts.

This patient representation allows us to convert from any common data model
(e.g. OMOP) to patient sequences and vice versa without any loss of temporal
information. To formulate this property, let’s denote $D_{i}$ to be data
associated with the $i$th patient in the source format, $P_{i}$ to be the
patient sequence converted from $D_{i}$, and $F$ to be the function that
converts source data to patient sequences, represented as $P_{i}=F(D_{i})$.
Let’s then denote $D_{i}^{\prime}$ to be the reconstructed source data and
$F^{\prime}$ to be the inverse function that converts patient sequences back
to original source format. This indicates $D_{i}^{\prime}=F^{\prime}(P_{i})$.
Finally, let’s denote $T$ to be a function that extracts all the dates for a
set of patient records. The patient representation is said to preserve
temporal information perfectly if and only if the following statement is true
for every single patient, $T(D_{i})=T(D_{i}^{\prime})+C_{i}$, where $C_{i}$ is
a constant that represents a consistent time shift e.g., $C_{i}=4$ days.

Figure 2: The patient representation preserves demographics, visit types, and
temporal intervals between visits and inpatient duration. It’s designed to
have the demographic prompt at the beginning including year at the first
visit, age at the first visit, gender and race tokens, then followed by a
series of visit blocks to represent the complete patient timeline. An
artificial time token (ATT) is inserted between the neighboring visit blocks
to keep track of the time intervals in days. In each visit block, all the
essential information is retained including the visit type and domain records.
In the case of inpatient visits, the inpatient ATT tokens (representing time
intervals in days) are inserted between groups of concepts that occur on the
same day, in addition, a discharge token is provided at the end of the visit
block.

###  3.2 OMOP Encoder

To create a patient sequence, we began by generating a demographic prompt
using data from the OMOP person and visit tables. This prompt included
essential demographic information such as the patient’s age at their initial
visit, the year of their first visit, their gender, and their race.
Subsequently, we constructed a series of visit blocks to represent the
patient’s entire medical history. We inserted an ATT token between these visit
blocks to signify the time intervals between them. Within each visit block, we
gathered all relevant records from OMOP domain tables (e.g., condition) and
arranged them chronologically based on their respective timestamps. In cases
where there are timestamp ties, we sorted the concepts. Additionally, three
artificial tokens (VS, VE, and VTT) were added at the beginning and end of
each visit block to denote the start, end, and type of the visit, as
illustrated in Figure 2. For inpatient visit blocks, extra processing steps
were necessary. Initially, we grouped records by their timestamps to identify
all inpatient spans, arranging them in chronological order. Next, we inserted
IATT between these spans, aligning them with the respective time intervals.
Finally, the discharge facility code was extracted from the OMOP visit table
and appended to the end of the block.

###  3.3 Generative Model

We used a GPT model with standard transformer decoders, where the input layer
utilized concept embedding and trainable positional embedding. The model was
trained using the Next Word Prediction learning objective. When generating a
patient sequence, we randomly sampled a demographic prompt from our source
sequences, which served as the input to the GPT model. Using these prompts,
the entire patient history was generated autoregressively by sampling tokens
from the predictive distribution at the final layer.

###  3.4 OMOP Decoder

The patient sequence was converted back to the OMOP format using the OMOP
decoder. The start-year prompt determined the EHR history’s beginning, using
January 1st as the default. Demographic data was stored in the person table,
while concepts were transformed into condition, drug, and procedure tables. A
date cursor was used to represent the “current time” as we were processing
each patient sequence, it was initially set to the star year and was updated
whenever an ATT token was encountered. We first parsed out the number of days
represented by the ATT token, and then moved the data cursor by the same
number of days to the future. During the sequence processing, the VS token
marked the start of a new visit block. We, therefore, extracted the token
corresponding to the visit type (the token immediately followed by VS) and
created a new visit record with the corresponding type.

All the tokens subsumed by this visit block were converted to condition, drug,
and procedure records and linked to the current visit. For outpatient visits,
we assumed that all data points were generated on the same day; therefore, we
set the start and end dates of the visit to the current value of the date
cursor. Similarly, all the domain records were set to the date cursor as well.

When processing inpatient visits, the date cursor was updated inside the visit
block based on the time intervals represented by IATT tokens between inpatient
spans. Domain records were generated with the current value of the date
cursor. Towards the end of the inpatient visit block, we extracted the last
token (the token right before VE) corresponding to a discharge facility and
updated the visit end date using the current value of the date cursor (as the
date cursor was frequently updated inside an inpatient visit block). This
allowed us to preserve the complete information about inpatient visits. If
generated sequences do not follow patterns presented in the patient
representation, they will be discarded to ensure the quality of the synthetic
data.

##  4 Experiments and Results

###  4.1 Data and Preprocessing

The source patient sequences were generated from the OMOP database derived
from Columbia University Irving Medical Center-New York Presbyterian Hospital
EHR data, which includes 3.7 million unique patients’ medical histories
including condition, medication, and procedure. Unknown concepts (i.e.,
$concept\\_id=0$) were removed from all domains except for the visit type when
constructing the patient sequences using the proposed patient representation.
Patients with less than 20 tokens were removed from the training dataset, and
approximately 2.3 million patients were included for training whereas 75,000
patients were held out for privacy evaluations. For the GPT model, we used a
context window of 512, 16 transformer decoders, 8 attention heads with a
dropout rate of 0.1, and 128 dimensions for both the embedding and hidden
units. All patients with longer than 512 tokens were post-truncated to fit the
context window. The statistics of training data was summarized in Table 1. We
trained the model for 2 epochs on 2 Nvidia 2080 TI GPUs with a batch size of
32 and a learning rate of 0.0002. The model checkpoint was created every
10,000 steps.

|  No. of visits per patient |  Sequence length per patient
---|---|---
mean |  16 |  148
std |  19 |  154
min |  2 |  20
25% |  4 |  38
50% |  8 |  78
75% |  21 |  198
max |  102 |  512

Table 1: Summary statistics of the CUIMC-NYP OMOP training data

During the first epoch, we used a standard data generator strategy, where
every training example was fed to the model; however, we switched to a random
sampling strategy to draw training examples during the second epoch. For
synthetic data generation, we used the 10th model snapshot because early
experiments showed its superior performance compared to other snapshots.
Furthermore, we used several sampling hyper-parameters including top k=300,
top k=200, top k=100, top p=90%, top=95%, and top=100% to generate different
synthetic OMOP datasets. Using specific top p/k values is a common technique
in language models for data generation. For instance, the top k approach
limits the selection of the k most probable tokens in the prediction
distribution during sampling. On the other hand, the top p method selects a
set of most likely tokens whose combined probability reaches p% (e.g. 90%) in
the predictive distribution. For each sampling strategy, 1M synthetic
sequences were generated and converted to OMOP. On average, 98% of the
generated sequences passed the validation and were converted to OMOP.

###  4.2 Evaluations

We followed the evaluation procedures proposed by [17] to compute the data
utility metrics including dimension-wise distribution, co-occurrence
relationship, and machine learning model performance. Some of the metrics were
originally designed for tabular EHR data; therefore, we adapted them to the
time-series setting. When using Kullback-Leibler (KL) divergence to evaluate
source and synthetic datasets, we use concept probabilities defined as,

$P_{prob}(c)=\frac{\sum_{i}^{n}\mathbbm{1}\Big{[}c\in
h_{i}\Big{]}}{\sum_{i}^{n}\sum_{j}^{m}\mathbbm{1}\Big{[}c_{j}\in
h_{i}\Big{]}}$

where $c$ denotes the target concept, $h_{i}$ denotes the $ith$ patient
history, $n$ and $m$ denote the total number of patients and concepts
respectively. Due to the small probability values, we opted to use the
prevalence instead for data visualization using a slightly modified formula
below,

$P_{prev}(c)=\frac{\sum_{i}^{n}\mathbbm{1}\Big{[}c\in h_{i}\Big{]}}{n}$

The only difference between the prevalence and the probability is the
normalization constant used as the denominator.

For baseline comparison, we added three variants of GPT models trained on
slightly different patient representations. The first baseline model was
trained on the adjusted CEHR-BERT representation, which differs from the
proposed patient representation by 1) CEHR-BERT ATT tokens contain a mix of
day/week/month/year tokens to represent time intervals while CEHR-GPT only
used the day tokens and 2) the IATT tokens and discharge facility tokens were
not used in the CEHR-BERT representation. The second baseline model used the
proposed patient representation with IATT tokens removed. We will refer to
this baseline as GPT-OUTPAT. The last baseline GPT was trained on the patient
representation widely used in time-series EHR research, where the concepts
were simply ordered chronologically and put in a sequence without any
additional artificial tokens. To use such sequences for comparison, we assumed
that all medical events in the patient sequence belonged to a single visit.
This baseline model will be referred to as GPT-Vanilla. We only used the
top_p=95% sampling strategy to generate synthetic OMOPs for these baseline
models. The comparison between these patient representations can be seen in
Supplementary Figure 9.

####  4.2.1 Dimension-wise Distribution

KL divergence was assessed to compare the concept probability distributions
between synthetic and source datasets among the entire population. In Figure
3, synthetic datasets generated by different patient representations and
sampling strategies were evaluated against real patient data. The results
showed that baseline models CEHR-BERT, GPT-Vanilla and GPT-OUTPAT with
sampling strategy using the threshold of top_p$=95\%$ diverged the least from
real concept probability distributions followed by GPT models trained with
top_p$=100\%$ and top_p$=95\%$. GPT models with threshold of top_k$=300$ and
top_k$=200$ with relatively similar divergence. However, GPT models with
threshold of top_k$=100$ sampling strategy had the largest KL divergence.

Figure 3: KL divergence for comparing concept probability distribution between
synthetic data and real data. The probabilities of concepts were calculated on
the scale of the entire population.

Furthermore, we conducted a qualitative analysis using the synthetic data
top_p=95% to gain a more comprehensive understanding. The dimension-wise
distributions between the synthetic and source datasets were compared at three
distinct levels: the entire population, specific sub-groups (e.g., female
population), and particular cohorts (e.g., hospitalization cases). Figure 4
illustrates the concept prevalence comparison between the original OMOP
dataset and the generated OMOP dataset, using a threshold of top_p$=95\%$. In
the high-frequency regions, most data points cluster closely around the
diagonal line, indicating a strong agreement between the source and synthetic
data. Conversely, data points appear more dispersed for low-frequency
concepts. Notably, in the subplot representing female conditions (located in
the first column, second row), there is an unusual cluster of concepts
positioned above the diagonal line. Further examination revealed that these
concepts were male-specific and should not appear in the female population.
Although there are a few instances of such cases in the source data, GPT
amplified such cases in synthetic data.

Figure 4: Concept prevalence comparison between the source OMOP and generated
OMOP using top $p=95\%$ in the log scale stratified by domain in columns and
by population in rows, where x-axis and y-axis represent the source and the
synthetic data respectively, and each dot represents a concept

In addition, we conducted a detailed comparison of the visit tables between
the two OMOP datasets. Our specific focus was on performing demographic
breakdowns and analyses for gender, race, and age group. The supplementary
Figures 10, 12, 11 highlight the top 10 most prevalent visits for each
demographic breakdown, showcasing that the trends in both the source and
generated visit tables exhibited notable similarities.

####  4.2.2 Co-occurrence Relationship

To measure how closely the generated datasets resemble the source, we computed
the KL divergence between their co-occurrence matrices. The matrix was
constructed temporally with the following logic: 1) for each concept in the
patient sequence, only the future concepts were taken into consideration for
creating concept pairs; 2) each patient could only contribute to the same
concept pair once; 3) the matrix was normalized into a proper probability
distribution by dividing the occurrences of each concept pair by the overall
number of pairs. Additionally, we set benchmarks for our analysis: a lower-
bound and an upper-bound. The lower bound was determined by applying the KL
divergence method to two random samples from the source data. The upper bound
was established by creating a theoretical co-occurrence matrix under the
assumption that all concepts in the source data were independent. The KL
divergence was then applied to this hypothetical matrix to calculate the upper
bound. Figure 5 illustrates that the datasets with CEHR-BERT and top_k=300
most closely approach the lower bound, with the top_p=95% and GPT-OUTPAT
baselines coming next. Datasets with top_k=100 and top_k=200 had marginally
higher KL divergence values, while the top_p=100%, the GPT-Vanilla datasets
exhibited the largest KL divergence.

To thoroughly examine the similarities between the original and the synthetic
data, we carried out a qualitative analysis using one of the synthetic
datasets. This involved comparing the most frequent pairs of co-occurring
concepts within each category, such as condition-condition (interpreted as one
condition concept followed by another condition).

Figure 5: KL divergence associated with different synthetic data. The closer
to the lower bound in the bottom left corner, the better the synthetic data.

Figure 6: Top 100 co-occurring concept pairs for each co-occurrence category
e.g. condition-condition interpreted as a condition concept followed by
another condition concept. The x and y axes represent the synthetic and source
data respectively

The results, illustrated in Figure 6, display the top 100 pairs for each type
of co-occurrence relationship. The analysis revealed that the synthetic data
accurately mirrored the co-occurrence patterns in most categories, with the
exception of the categories ending with a drug concept (as shown in the second
column of Figure 6), where the data points were more scattered. Notably, a
majority of the points in condition-condition and procedure-procedure pairs
aligned along the diagonal line.

Figure 7: Comparison of temporal co-occurrence networks around Type 2 Diabetes
Mellitus (T2DM) for real (left) and synthetic data (right). Networks were
constructed with T2DM’s (large circle) top 5 co-occurring condition concepts
and their own top 5 co-occurring condition concepts. The metrics Prevalence
(top) and Pointwise Mutual Information-3 (PMI3) (bottom) were used to quantify
co-occurrence. Gold edges indicate edges that are shared by both real and
synthetic data, edge thickness indicates the strength of co-occurrence, and
arrow direction indicates the direction of the temporal association.

For an in-depth analysis, we focused on examining the co-occurrence
relationships associated with Type 2 Diabetes Mellitus (T2DM). Figure 7
qualitatively compares the real and top_p=95% synthetic data networks of
condition concepts that co-occur around T2DM. The degree of co-occurrence is
calculated with prevalence, which evaluates the frequency of concept co-
occurrence, and Pointwise Mutual Information-3 (PMI3), which evaluates the
probabilistic association of concepts. There is much overlap between the real
and synthetic networks when prevalence is used to construct them, but
connections are more disparate when PMI3 is used.

####  4.2.3 Predictive performance

For this analysis, we constructed five prediction tasks using the method
described in [13] and the Book of OHDSI [18]. Table 2 shows the cohort and the
corresponding definition.

Cohort |  Definition
---|---
HF readmission |  HF patients who have a 30-day all-cause readmission. Observation window: 360 days, Prediction windows 30 days
Hospitalization |  2-year risk of hospitalization starting from the 3rd year since the initial entry into the EHR system. Observation window: 540 days, hold-off window: 180 days, prediction windows: 720 days
COPD readmission |  COPD patients who have a 30-day all-cause readmission. Observation window: 360 days, prediction windows: 30 days
Afib ischemic stroke |  Afib patients with 1-year risk since the initial diagnosis of afib ischemic stroke. Observation window: 720 days, prediction windows: 360 days
CAD CABG |  Patients initially diagnosed with Coronary Arterial Disease (CAD) without any prior stent graft will receive the Coronary artery bypass surgery (CABG) treatment. Observation window: 720 days, prediction windows: 360 days
Table 2: Cohort definitions

To extract features, we first rolled up the medical concepts using ontological
hierarchies to reduce dimensionality [see supplementary materials] and used
the bag-of-word (BOW) approach, where we counted the frequency of each concept
in a given observation window. For each task, we split the cohort data (both
synthetic and real) into training and testing sets with a split ratio of
85:15. We ran logistic regression using Sklearn’s implementation with the
default configuration. Finally, the area under the receiver operating
characteristics curve (AUC) was calculated using the test set. In addition, we
reported PR-AUC (precision-recall) due to the class imbalance often present in
EHR data. Table 3 shows the prevalence of the positive cases, ROC-AUC, and PR-
AUC for each synthetic data. The metrics associated with the baseline CEHR-
BERT are reported in the Supplementary Materials Table 7. For easy comparison
of synthetic datasets, we defined a consolidated distance metric as the
weighted average of the relative differences of the three aforementioned
metrics,

$dist=\frac{|\delta_{Pre}|\times 0.5}{Pre_{true}}+\frac{|\delta_{AUC}|\times
0.25}{AUC_{true}}+\frac{|\delta_{PR}|\times 0.25}{PR_{true}}$

where $\delta_{Pre}$, $\delta_{AUC}$, and $\delta_{PR}$ represent the
differences in prevalence, ROC-AUC, and PR-AUC between the source and the
synthetic data; $Pre_{true}$, $AUC_{true}$, and $PR_{true}$ denote the ground
truth metrics generated from the source data.

Figure 8 presents the distance metrics for various synthetic datasets across
different cohorts. The dataset created with top_k=300 displayed the best
performance in HF and COPD Readmission but was less effective in other
cohorts. The top_p=95% dataset maintained consistent performance levels,
ranking as the best in Hospitalization and Afib ischemic stroke, and as
second-best in HF and COPD readmission. In comparison, the GPT-OUTPAT baseline
dataset exhibited a slightly higher divergence across the Hospitalization,
Afib Ischemic Stroke, and CAD CABG groups. However, the CEHR-BERT baseline
showed similar patterns in most cohorts, with a notably lower divergence in
CAD CABG. It is noted that both GPT-OUTPAT and CEHR-BERT were not included in
the HF Readmission and COPD Readmission analyses due to the absence of
synthetic patients meeting the selection criteria for these cohorts. Other
synthetic datasets, specifically those generated with top_p=100% and top_k=200
lagged behind in performance. Interestingly, top_k=100 showed a unique
pattern, being the closest in the distance for the CAD CABG cohort but
underperforming in the others.

Figure 8: The consolidated distance metrics for different synthetic datasets
stratified by cohort. GPT-Vanilla and CEHR-BERT were omitted as the cohorts
couldn’t be constructed due to the loss of temporal information. GPT-OUTPAT
was omitted from HF readmission and COPD readmission as the cohorts for the
same reason.

Cohort | Real | p=95% | p=100% | k=100 | k=200 | k=300 | GPT-OUTPAT
---|---|---|---|---|---|---|---
HF readmission |  Pre = 25.7 AUC = 65.7 PR = 39.3 |  Pre = 27.6 AUC = 69.2 PR = 45.7 |  Pre =27.7 AUC = 52.4 PR = 29.0 |  Pre = 30.7 AUC = 68.1 PR = 47.8 |  Pre = 29.3 AUC = 54.0 PR = 32.9 |  Pre = 26.5 AUC = 61.1 PR = 33.8 |  Pre =100.0 AUC = NA PR = NA
Hospitalization |  Pre = 5.6 AUC = 75.3 PR = 19.5 |  Pre = 5.2 AUC = 77.1 PR = 21.4 |  Pre = 7.4 AUC = 71.3 PR = 20.2 |  Pre = 2.8 AUC = 87.0 PR = 22.1 |  Pre = 5.2 AUC = 84.2 PR = 20.8 |  Pre = 6.3 AUC = 78.7 PR = 24.6 |  Pre = 5.2 AUC = 70.2 PR = 14.3
COPD readmission |  Pre = 34.5 AUC = 74.2 PR = 83.8 |  Pre = 37.8 AUC = 76.4 PR = 84.4 |  Pre = 47.2 AUC = 74.1 PR = 67.2 |  Pre = 26.4 AUC = 75.9 PR = 90.3 |  Pre = 28.3 AUC = 70.1 PR = 82.8 |  Pre = 34.5 AUC = 68.8 PR = 80.2 |  Pre = NA AUC = NA PR = NA
Afib ischemic stroke |  Pre = 8.7 AUC = 84.0 PR = 48.5 |  Pre = 10.2 AUC = 78.9 PR = 41.2 |  Pre = 10.4 AUC = 70.7 PR = 39.1 |  Pre = 16.6 AUC = 77.1 PR = 50.5 |  Pre = 15.8 AUC =68.9 PR = 36.6 |  Pre = 10.8 AUC = 76.8 PR = 38.5 |  Pre = 9.7 AUC = 67.2 PR = 27.2
CAD CABG |  Pre = 7.1 AUC = 88.4 PR = 55.9 |  Pre = 4.1 AUC = 81.5 PR = 25.2 |  Pre = 4.4 AUC = 52.9 PR = 4.3 |  Pre = 7.2 AUC = 84.7 PR = 31.3 |  Pre = 4.9 AUC = 73.5 PR = 24.3 |  Pre = 4.0 AUC = 79.0 PR = 24.1 |  Pre = 3.5 AUC = 81.5 PR = 44.4
Table 3: Logistic regression model performance across different datasets. In
each cell, three numbers were reported including the prevalence of the
positive cases, ROC-AUC, and PR-AUC

###  4.3 Privacy Evaluations

We adopted the privacy evaluation framework outlined in [17], focusing on
quantifying the risk of privacy breaches through attribute and member
inference attacks. The attribute inference risk assesses the likelihood of
inferring sensitive attributes from a synthetic dataset by utilizing targeted
patients’ demographic information and common diagnoses found in real datasets
available to the adversary. In a membership inference attack, the adversary
tries to determine if the target record was in the training set given all or
partial attributes.

####  4.3.1 Membership Inference Attack

We simulated two different attack scenarios, namely the dataset attack and the
model attack. In a dataset attack, attackers only have access to the synthetic
dataset, whereas they can query the model itself in a model attack. Their goal
is to find a way to infer whether or not a real patient record was used for
training [19]. For example, the model might output a high probability
associated with the training data points and a low probability for the non-
training ones, similarly, the synthetic data might be more similar to the
training data than the non-training ones.

The attack dataset was constructed based on the following procedure 1) we
created the negative examples by taking the 75000 holdout set and labeling
them as negative; 2) we created the positive examples by randomly selecting
75000 patients from the training set and assigning them as positive. In the
dataset attack, for each attack data point, we found the best match from the
synthetic data using a hamming distance metric, which sums up the absolute
differences in age, number of visits, and concepts [19], finally, the median
distance was used as a threshold, below which patients were predicted as
positive and negative otherwise. In the model attack, we computed the average
loss of each patient sequence by feeding the attack dataset to the model
directly, similarly, the median loss was used to infer the positive/negative
cases, where the positive cases were assigned to those below the threshold and
positive otherwise. Finally, we calculated recall/precision/f1/accuracy using
the ground truth labels against those inferred predictions. Table 4 shows that
the accuracy of both attacks is slightly less than 50%, indicating that the
performance of such attacks is worse than a random guess.

| Accuracy | Recall | Precision | F1
---|---|---|---|---
Model Attack | 0.4955 | 0.4955 | 0.4942 | 0.4949
Data Attack | 0.4941 | 0.4996 | 0.4959 | 0.4978

Table 4: Membership Inference Attack metrics

####  4.3.2 Attribute Inference Attack

In attribute inference attack, we assumed the attacker has access to a group
of target patients along with their demographic data, including age, gender,
race, the year of their initial clinical visit, and prevalent clinical
conditions like hypertension, abdominal pain, chest pain, etc. With these
features, the attacker identifies the patient in the synthetic dataset with
the closest attribute resemblance to the target individual. Then the attacker
aims to use the sensitive attributes of matched patient from synthetic dataset
to infer the corresponding sensitive attributes of the target patient.

To quantitatively evaluate the attribute inference risk, we randomly sampled
150,000 patients from the real patient dataset as a target group. The search
group compromised 1 million generated synthetic patients. To find most similar
patients from the search group, we created a set of common attributes
including the demographic data along with the top 1% of the most prevalent
condition concepts, represented as one-hot encoded features. Then a k-nearest
neighbors (KNN with k = 1) algorithm was applied to each target patient and a
synthetic patient with the smallest euclidean distance was found. Finally we
extracted the sensitive attributes (condition concepts not in the top 1% tier)
from matched target and synthetic patients. F1 scores were computed for the
sensitive attributes of each matched patient pair and aggregated across all
matched patients.

A baseline analysis was performed substituting real patients for synthetic
patients. A result lower than the baseline suggests that the likelihood of
finding a synthetic data similar to a real patient is lower than finding a
real patient similar to a real patient who share sensitive attributes. This
implies a lower attribute inference risk which could be acceptable. Our
training set was randomly divided into two halves, with 1 million real
patients assigned as the target group, and the remaining 1 million as the
search group. The matching process was consistently applied and an aggregated
F1 score was computed. Table 5 shows that the F1 score of synthetic vs real is
less than real vs real scenario. But it’s still higher than other models who
didn’t capture the real patterns as effectively. [19]

| Recall | Precision | F1
---|---|---|---
Synthetic vs Real | 0.0350 | 0.0343 | 0.0271
Real vs Real | 0.0612 | 0.0468 | 0.0421

Table 5: Attribute Inference Attack metrics

##  5 Discussion

To the best of our knowledge, this is the first attempt to utilize GPT for
generating time-series heterogeneous EHR data, while preserving patient
privacy. Our main contribution lies in designing a novel patient
representation that preserves a complete timeline of the patient’s history,
along with crucial visit details, thereby enabling GPT to create realistic
patient sequences. Importantly, this representation facilitates seamless
conversion back to the OMOP format, simplifying dissemination and analysis.
This patient representation could serve as an effective messenger for
transferring information across various standard data models. At present, our
system is tailored to the OMOP format. However, it is designed with
adaptability in mind, enabling us to seamlessly integrate new encoder/decoder
pairs. This flexibility would facilitate the conversion of patient sequences
to other widely-used data models, such as i2b2[20].

The study undertakes a three-tiered evaluation approach, systematically
comparing synthetic and real datasets based on their marginal (column-wise
distribution), conditional (co-occurrence relationship), and joint
distributions (predictive performance). Concurrently, as the evaluation
progresses through these three levels, there is a corresponding escalation in
the complexity and challenge of the tasks involved. The outcomes of the KL
divergence analysis revealed a nuanced relationship between the top k/p
sampling strategies and the performance across the evaluation levels.
Specifically, an increase in top k/p values enhanced performance in level 1
concept prevalence. However, an excessively high or low top k/p value
adversely affects both level 2 co-occurrence metrics and level 3 machine
learning predictions. This pattern suggests that including more tokens in the
predictive distribution introduces greater uncertainty and a wider array of
potential patient trajectory variations in the data generation process,
thereby escalating the difficulty of achieving comparable performance
outcomes.

The sampling strategies of top_p=95% and top_k=300 seem to be most effective
for generating the synthetic data. For instance, the synthetic data created
with top_p=95% demonstrates the second smallest divergence in both dimension-
wise distribution and co-occurrence relationship. Simultaneously, the
corresponding synthetic cohorts successfully replicated the performance
metrics in all predictive tasks, with the exception of CAD CABG. Finally, the
significance of this patient representation transcends synthetic data
generation; we believe it has the potential to establish the groundwork for
integrating time into patient representations across diverse EHR-based deep-
learning models.

###  5.1 Loss of Temporal Information

The reason that CEHR-GPT replicated the performance metrics of the prediction
tasks can be attributed to the use of time tokens in its underlying patient
representation. The majority of the prediction problems are phrased as “For a
group of target patients who share similar characteristics, who would
experience a particular medical event in one year from the index event?” in
EHR research, therefore maintaining a complete patient timeline is crucial for
time-sensitive cohort constructions [18].

We claimed that the proposed patient representation had almost zero loss of
temporal information, although this makes intuitive sense, there does not
exist a formal metric to quantify this. To bridge this gap, we conceived a new
metric named loss of temporal information (LOTI) to estimate the shrinkage of
the patient timeline due to the use of the patient representation in an EHR
dataset. Let’s denote $T$ to be the time interval measured in days, ATT to be
an artificial time token that represents a time interval ($W_{0}$), $F$ to be
a function that maps $T$ to an ATT token (four days to $W_{0}$). In addition,
let $G$ be the inverse function of $F$ that converts an ATT to $T$, moreover,
we impose the constraint on G such that it takes the lower bound of ATT e.g.,
$W_{0}\implies 0$ days. Formally, we define LOTI as the expected difference
between the original time interval $T$ and the reconstructed time interval
$G(F(T))$ as the following,

$LOTI=E_{p(T)}\Big{[}T-G\big{(}F(T)\big{)}\Big{]}$

where $P(T)$ is the probability of $T$ observed in the training data defined
as,

$P(T)=\frac{\textit{Freq of T}}{\Sigma\textsuperscript{T}\textit{Freq of T}}$

We computed LOTI for the patient representations utilized in CEHR-GPT and
baseline models, shown in Figure 9. As Table 6 shows, CEHR-GPT has the least
LOTI compared to the other patient representations while GPT-OUTPAT has a
slightly higher time shrinkage because the inpatient duration was not
retained. CEHR-BERT has a relatively large LOTI compared to the previous two
representations due to the use of coarse ATT tokens. Finally, GPT-Vanilla has
the most LOTI, which is equal to the expected length of the timeline in the
source population due to the complete collapse of the timeline.

Representation | Between visit ATT token | Between inpatient span ATT token | LOTI
---|---|---|---
CEHR-GPT |  Day token for $T\leq 1080$ LT token for $T>1080$ | Day token | 7.739
GPT-OUTPAT |  Day token for $T\leq 1080$ LT token for $T>1080$ | N/A | 7.962
CEHR-BERT |  Day token for $T<7$ Week token for $7\leq T<30$ Month token for $30\leq T<360$ LT token $T\geq 360$ | N/A | 31.482
GPT-Vanilla | N/A | N/A | 111.164
Table 6: Loss of Temporal Information for Different Patient Representations

###  5.2 Time Invariance and Sensitivity

Across all synthetic datasets, the dimension-wise distribution (marginal
distribution) and co-occurrence relationship were well preserved, regardless
of the patient timeline’s integrity within the models used. Even with a high
LOTI, baseline models such as CEHR-BERT accurately mirrored both marginal
distribution and co-occurrence relationship, yielding results similar to those
from CEHR-GPT. This implies that these two measures may be largely time-
invariant, unaffected by any shrinkage in the patient timeline. The rationale
behind this is rooted in their construction methods, which either disregard or
marginalize the temporal factor. Marginal distribution was constructed by
counting the unique number of patients associated with the target concept,
which was then normalized by a constant. The construction disregarded
temporality, as the placement of a concept on the timeline was not a factor of
consideration. The co-occurrence matrix was created using time initially but
was marginalized after all co-occurring pairs were collected from the patient
population.

On the contrary, the predictive performance is extremely sensitive to any
change made to the patient timeline as the cohort construction requires the
integrity of the patient timeline. For example, HF readmission requires a
30-day prediction window from the index event (defined as the hospitalization
episode with a heart failure diagnosis). Any shrinkage to the patient timeline
will disrupt the construction of this cohort. The synthetic HF readmission
cohort produced by CEHR-BERT showed a readmission rate of 100% due to the
shrinkage of the timeline in this patient representation; whereas, the actual
expected rate of readmission should be approximately 25%. Compared to CEHR-
GPT, GPT-OUTPAT encoded time intervals between visits but did not preserve the
duration of inpatient visits, therefore having a slightly higher LOTI. As a
consequence, it showed reasonable performance in the cohorts (Hospitalization
and Afib ischemic stroke), which had a 360-day prediction window and were thus
less impacted by timeline shrinkage. However, in the case of HF Readmission
(where a 100% readmission rate was observed) and COPD Readmission (which
identified no patients), GPT-OUTPAT was less successful. These cohorts used a
short 30-day prediction window, making them highly sensitive to any
distortions in the timeline, which likely led to synthetic patients not
meeting the cohort selection criteria. Interestingly, the CAD CABG cohort
presents a notable deviation from the general trend, where the CEHR-GPT
dataset with top_k=100 outperformed both the top_p=95% and top_k=300
configurations. Additionally, the CEHR-BERT synthetic data accurately
replicated the machine learning performance metrics as well. This indicates
that the CAD CABG cohort was less affected by time shrinkage.

Therefore, selecting the appropriate patient representation is pivotal in
maintaining specific properties of the source data when generating synthetic
data. The choice hinges on the intended application of synthetic data,
ensuring that critical features and patterns inherent to the patient
information are accurately reflected and retained.

###  5.3 Time Sensitive Forecasting

Because the patient representation encodes all the temporal information in the
sequence, the trained GPT model could be used potentially for time-sensitive
forecasting. We could prompt the trained GPT model with a patient history and
estimate the time of the next visit via a Monte Carlo Sampling approach shown
in the following equation,

$P(\delta_{t}|h)\approx\frac{\sum^{n}_{i=1}\mathbbm{1}\Big{[}M_{gpt}(h)=\delta_{t}\Big{]}}{n}$

where $M_{gpt}$ denotes the GPT model, $h$ denotes a patient history,
$\delta_{t}$ denotes any time interval, and $n$ represents the number of
samples. Then we can use the expectation $E\big{(}\delta_{t}\big{)}$ as the
predicted time interval. In addition, we can quantify the confidence by
calculating the standard deviation e.g. $sd(\delta_{t})$. Similarly, we can
predict the visit type ($v$) using the same Monte Carlo approach.

$P(v|E\big{(}\delta_{t}\big{)},h)\approx\frac{\sum^{n}_{i=1}\mathbbm{1}\Big{[}M_{gpt}\Big{(}E\big{(}\delta_{t}\big{)},h\Big{)}=v\Big{]}}{n}$

Finally, we can predict the most likely medical event ($c$) given the
predicted visit type $v$,

$P(c|v,E\big{(}\delta_{t}\big{)},h)\approx\frac{\sum^{n}_{i=1}\mathbbm{1}\Big{[}M_{gpt}\Big{(}v,E\big{(}\delta_{t}\big{)},h\Big{)}=c\Big{]}}{n}$

This approach goes beyond conventional prediction methods by not only
forecasting future medical events but also determining the timing of the next
visit and the specific medical events associated with that visit type. This
predictive model could provide a more detailed and actionable timeline for
patient care.

###  5.4 Limitations

While synthetic datasets demonstrated a high degree of similarity to source
data, they are subject to certain known constraints. Firstly, a selection bias
was present in the training data due to the constraints on sequence length,
ranging from 20 to 512. This limitation resulted in the partial inclusion of
patients with chronic conditions, which typically require a longer context
window for accurate representation. While extending the context window of the
model could potentially address this issue, it may introduce unforeseen
effects. Finding the optimal configuration to accommodate a broader context
window would require comprehensive experiments.

Secondly, identifying an optimal sampling strategy for generating synthetic
data remains a challenge due to the presence of numerous hyperparameters such
as temperature, top_p, and top_k. These parameters, when used in conjunction,
could yield a wide array of configurations. While the top_p=95% strategy
showed the least divergence from the source data, it was unable to accurately
replicate performance metrics for CAD CABG. As an interim solution, we may
publish multiple versions of synthetic data, along with their corresponding
performance metrics. This approach would allow researchers to select the most
suitable dataset for their specific use case.

Thirdly, the GPT model showed a propensity to over-represent prevalent
concepts, skewing towards those with higher frequencies in the dataset. An
illustration of this is seen in the synthetic data, where 78% of patients had
at least one outpatient visit, in contrast to the actual data where this
figure was 73%. This discrepancy indicates a bias in the model towards more
common occurrences. In addition, the co-occurrence analysis also demonstrated
that the synthetic reconstruction faithfully represents the frequent concept
pairs in the original data but may be less effective at recovering the finer
associations between rare concepts as shown in Figure 7. To address the over-
representation of prevalent concepts by the GPT model, future work will look
into implementing regularization techniques. One promising approach is
adaptive regularization, which can be outlined in several steps: 1) Model
Training: Begin by training the GPT model for a predetermined number of steps;
2) Sequence Generation and Analysis: Generate a sample of patient sequences
from the trained model and calculate the marginal distribution of the concepts
within these sequences; 3) Distribution Comparison and Adjustment Score
Calculation: Compare the model-generated distribution to the empirical
distribution derived from the actual data. From this comparison, calculate an
adjustment score for each concept; 4) Logit Adjustment: Modify the logits for
each concept in the model according to the calculated adjustment scores. By
implementing this procedure, the influence of each concept on the model’s
learning process would be adaptively modified during back-propagation,
allowing for an update in the model parameters that takes into account the
disparity between the generated and actual data distributions. This should
help in reducing the bias towards over-represented concepts.

Furthermore, improving patient representation is an area for further
development. Currently, the model’s representation is limited to daily
intervals and does not capture more precise measurements like hours or
minutes. This limitation is particularly relevant for intensive care unit
(ICU) data where time-sensitive decisions are critical. Furthermore, the
current model framework assigns the first visit of every synthetic patient to
the start of the year, as it only includes a year token to denote the
commencement of patient history. To refine the accuracy of patient history
commencement, integrating a month token alongside the year token is being
considered to properly represent seasonality. This would provide a more
accurate reconstruction of the starting point for a patient’s first visit in
the generated data.

Lastly, there is a necessity to incorporate the death event within patient
sequences. Including this event would allow the synthetic data to more
accurately represent mortality, enhancing its utility for predictions related
to patient outcomes and lifespan. These enhancements aim to create a more
precise and clinically relevant synthetic dataset that better mirrors the
complexities of real-world patient trajectories.

##  6 Conclusion

To our knowledge, this is the first attempt to utilize GPT for generating
time-series EHR data. Our main contribution lies in the design of a patient
representation that captures temporal dependencies among token types, enabling
GPT to generate realistic patient sequences. Moreover, this representation
allows for easy conversion back to the OMOP format. Comprehensive evaluations
showed that the synthetic data effectively captures the intricate patterns
present in EHR data.

## References

  * [1] Murray, R. E., Ryan, P. B. & Reisinger, S. J.  Design and validation of a data simulation model for longitudinal healthcare data.  In _AMIA Annual Symposium Proceedings_ , vol. 2011, 1176 (American Medical Informatics Association, 2011).
  * [2] Ghosheh, G., Li, J. & Zhu, T.  A review of generative adversarial networks for electronic health records: applications, evaluation measures and data sources (2022).  URL https://arxiv.org/abs/2203.07018v2.
  * [3] Choi, E. _et al._ Generating multi-label discrete patient records using generative adversarial networks.  In Doshi-Velez, F. _et al._ (eds.) _Proceedings of the 2nd Machine Learning for Healthcare Conference_ , vol. 68 of _Proceedings of Machine Learning Research_ , 286–305 (PMLR, 2017).  URL https://proceedings.mlr.press/v68/choi17a.html.
  * [4] Li, Y., Bengio, S. & Du, N.  Time-dependent representation for neural event sequence prediction.  _6th International Conference on Learning Representations, ICLR 2018 - Workshop Track Proceedings_ (2017).  URL https://arxiv.org/abs/1708.00065v4.
  * [5] Li, J., Cairns, B. J., Li, J. & Zhu, T.  Generating synthetic mixed-type longitudinal electronic health records for artificial intelligent applications.  _npj Digital Medicine 2023 6:1_ 6, 1–18 (2023).  URL https://www.nature.com/articles/s41746-023-00834-7.
  * [6] Lee, D. _et al._ Generating sequential electronic health records using dual adversarial autoencoder.  _Journal of the American Medical Informatics Association : JAMIA_ 27, 1411–1419 (2020).  URL https://pubmed.ncbi.nlm.nih.gov/32989459/.
  * [7] Hripcsak, G. & Albers, D. J.  Next-generation phenotyping of electronic health records.  _Journal of the American Medical Informatics Association : JAMIA_ 20, 117–121 (2013).  URL https://pubmed.ncbi.nlm.nih.gov/22955496/.
  * [8] Openai, A. R., Openai, K. N., Openai, T. S. & Openai, I. S.  Improving language understanding by generative pre-training URL https://gluebenchmark.com/leaderboard.
  * [9] Vaswani, A. _et al._ Attention is all you need .
  * [10] Goodfellow, I. _et al._ Generative adversarial networks.  _Communications of the ACM_ 63, 139–144 (2014).  URL https://arxiv.org/abs/1406.2661v1.
  * [11] Zhang, Z., Yan, C., Lasko, T. A., Sun, J. & Malin, B. A.  Synteg: a framework for temporal structured electronic health data simulation.  _Journal of the American Medical Informatics Association_ 28, 596–604 (2021).  URL https://dx.doi.org/10.1093/jamia/ocaa262.
  * [12] Gulrajani, I., Ahmed, F., Arjovsky, M., Dumoulin, V. & Courville, A.  Improved training of wasserstein gans.  _Advances in Neural Information Processing Systems_ 2017-Decem, 5768–5778 (2017).  URL https://arxiv.org/abs/1704.00028v3.
  * [13] Pang, C. _et al._ Cehr-bert: Incorporating temporal information from structured ehr data to improve prediction tasks.  In Roy, S. _et al._ (eds.) _Proceedings of Machine Learning for Health_ , vol. 158 of _Proceedings of Machine Learning Research_ , 239–260 (PMLR, 2021).  URL https://proceedings.mlr.press/v158/pang21a.html.
  * [14] Li, Y. _et al._ Behrt: Transformer for electronic health records.  _Scientific Reports 2020 10:1_ 10, 1–12 (2020).  URL https://www.nature.com/articles/s41598-020-62922-y.
  * [15] Rasmy, L., Xiang, Y., Xie, Z., Tao, C. & Zhi, D.  Med-bert: pretrained contextualized embeddings on large-scale structured electronic health records for disease prediction.  _npj Digital Medicine 2021 4:1_ 4, 1–13 (2021).  URL https://www.nature.com/articles/s41746-021-00455-y.
  * [16] Kraljevic, Z. _et al._ Foresight – generative pretrained transformer (gpt) for modelling of patient timelines using ehrs (2022).  URL https://arxiv.org/abs/2212.08072v2.
  * [17] Yan, C. _et al._ A multifaceted benchmarking of synthetic electronic health record generation models.  _Nature Communications_ 13 (2022).  URL https://doi.org/10.1038%2Fs41467-022-35295-1.
  * [18] OHDSI.  _The Book of OHDSI: Observational Health Data Sciences and Informatics_ (OHDSI, 2019).  URL https://books.google.com/books?id=JxpnzQEACAAJ.
  * [19] Theodorou, B., Xiao, C. & Sun, J.  Synthesize high-dimensional longitudinal electronic health records via hierarchical autoregressive language model.  _Nature Communications 2023 14:1_ 14, 1–13 (2023).  URL https://www.nature.com/articles/s41467-023-41093-0.
  * [20] Murphy, S. N. _et al._ Serving the enterprise and beyond with informatics for integrating biology and the bedside (i2b2).  _Journal of the American Medical Informatics Association_ 17, 124–130 (2010).

##  7 Supplementary Materials

Figure 9: The comparison of the patient representations.

Cohort | Real | GPT-OUTPAT | CEHR-BERT
---|---|---|---
HF readmission |  Pre = 25.7 AUC = 65.7 PR = 39.3 |  Pre =100.0 AUC = NA PR = NA |  Pre =100.0 AUC = NA PR = NA
Hospitalization |  Pre = 5.6 AUC = 75.3 PR = 19.5 |  Pre = 5.3 AUC = 70.2 PR = 14.3 |  Pre = 0.3 AUC = 76.0 PR = 1.0
COPD readmission |  Pre = 34.5 AUC = 74.2 PR = 83.8 |  Pre = NA AUC = NA PR = NA |  Pre = NA AUC = NA PR = NA
Afib ischemic stroke |  Pre = 8.7 AUC = 84.0 PR = 48.5 |  Pre = 9.7 AUC = 67.2 PR = 27.2 |  Pre = 19.5 AUC = 72.4 PR = 60.0
CAD CABG |  Pre = 7.1 AUC = 88.4 PR = 55.9 |  Pre = 3.5 AUC = 81.5 PR = 44.4 |  Pre = 5.4 AUC = 77.2 PR = 44.1
Table 7: Logistic regression (LR) model performance across baseline synthetic
datasets. In each cell, three numbers were reported including the prevalence
of the positive cases, ROC-AUC, and PR-AUC. If the prevalence=NA, this
indicates that $0$ patients were identified in the cohort. AUC=NA and PR=NA
indicate the LR model could not be run successfully due to either 100% or 0%
prevalence.

Figure 10: The visit prevalence stratified by age group

Figure 11: The visit prevalence stratified by race

Figure 12: The visit prevalence stratified by gender


Ronald A. Remmerswaal

Arthur E.P. Veldman

Bernoulli Institute, University of Groningen

PO Box 407, 9700 AK Groningen, The Netherlands

two-phase flow volume of fluid method parabolic reconstruction

This work is part of the research programme SLING, which is (partly) financed by the Netherlands Organisation for Scientific Research (NWO).


# Formation of moons and equatorial ridge around top-shaped asteroids after
surface landslide

Ryuki Hyodo ISAS/JAXA, Sagamihara, Kanagawa, Japan<EMAIL_ADDRESS>Keisuke Sugiura Earth-Life Science Institute, Tokyo Institute of Technology,
Meguro-ku, Tokyo 152-8550, Japan

###### Abstract

Top-shaped asteroids have been observed among near-Earth asteroids. About half
of them are reported to have moons (on the order of $\sim 1$wt.% of the top-
shaped primary) and many of them have an equatorial ridge. A recent study has
shown that the enigmatic top-shaped figure of asteroids (e.g., Ryugu, Bennu,
and Didymos) could result from an axisymmetric landslide of the primary during
a fast spin-up near the breakup rotation period. Such a landslide would
inevitably form a particulate disk around an asteroid with a short timescale
($\sim 3$ hours). However, the long-term full dynamical evolution is not
investigated. Here, we perform a continuous simulation ($\sim 700$ hours) that
investigates the sequence of events from the surface landslide that forms a
top-shaped asteroid and a particulate disk to disk evolution. We show that the
disk quickly spreads and produces moons (within $\sim 300$ hours). The mass of
the formed moon is consistent with what is observed around the top-shaped
asteroids. We also demonstrate that an equatorial ridge is naturally formed
because a fraction of the disk particles re-accretes selectively onto the
equatorial region of the primary. We envision that Ryugu and Bennu could once
have an ancient moon that was later lost due to a successive moon’s orbital
evolution. Alternatively, at the top-shaped asteroid that has a moon, such as
Didymos, no significant orbital evolution of the moon has occurred that would
result in its loss. Our study would also be qualitatively applicable to any
rubble-pile asteroids near the breakup rotation period.

Asteroids (72); Near-Earth objects (1092)

##  1 Introduction

Figure 1: Schematic summary of our paper. Panel (a): a rubble-pile asteroid
spins up due to, for example, the YORP effect, small impacts, a close
encounter with a planet, or re-accumulation after a catastrophic impact. Panel
(b): a surface landslide occurs when a critical spin-state is realized and a
top-shaped figure is formed (Sugiura et al., 2021). Panel (c): a particulate
disk spreads due to inelastic collisions and gravitational interactions among
particles. Panel (d): a moon is gravitationally accreted outside the Roche
limit of the central top-shaped body and an axisymmetric equatorial ridge is
formed due to the re-accretion of disk particles. Panels (e) and (f): the
formed moon is lost or remains, depending on the long-term orbital evolution
between the moon and the primary.

More and more growing interests are on near-Earth asteroids (NEAs) together
with recent activities of asteroid exploration missions. Top-shaped asteroids
may exist ubiquitously among NEAs (Walsh & Jacobson, 2015; Margot et al.,
2015). These include asteroid 162173 Ryugu (Watanabe et al., 2019), 101955
Bennu (Lauretta et al., 2019), 65803 Didymos, 1999 KW4 (Ostro et al., 2006),
and 2001 SN263 (Becker et al., 2015). Their diameters are less than several
kilometers, and recent in-situ observations by spacecrafts showed that they
might be rubble-pile bodies (e.g., Watanabe et al., 2019). Because of their
small size, Yarkovsky-O’Keefe-Radzievskii-Paddack (YORP) effect (Rubincam,
2000), small impacts (Takeda & Ohtsuki, 2009), and/or a close encounter with a
planet (Hyodo et al., 2016), for example, could efficiently change their spin
states. The acceleration rate of spin-up (or spin-down) depends on each
physical process (e.g., roughly instantaneous by an impact and $\gtrsim 100$
Kyr by YORP effect on a kilometer-sized body).

Interestingly, about half of the reported top-shaped asteroids have moons
around them (e.g., 2001 NE263; Becker et al. 2015, 1999 KW4; Ostro et al.
2006, 1994 CC; Brozović et al. 2011, Didymos; Naidu et al. 2020). The masses
of the moons are roughly on the order of $\sim 1$wt.% of the host top-shaped
asteroids. Furthermore, many of the top-shaped asteroids are reported to have
an equatorial ridge (Benner et al., 2015).

There are several proposed mechanisms to form a top-shaped asteroid. The top-
shaped figure may be formed during re-accretion after catastrophic disruption
of a parent body (Michel et al., 2020). Alternatively, a mass movement with
reshaping or a landslide of surface materials of a rubble-pile asteroid due to
fast spin-up may form the top-shaped figure (Walsh et al., 2008, 2012; Harris
et al., 2009; Hirabayashi et al., 2020; Sugiura et al., 2021).

Recently, Sugiura et al. (2021) used the Smoothed Particle Hydrodynamics (SPH)
method to study rotational deformation during spin-up, which has also been
used to model the shape deformation due to impacts (Jutzi & Asphaug, 2015). In
their study, they included effective bulk friction of rubble pile bodies as a
parameter (i.e., effective friction angle $\phi_{\rm fri}$), which could
effectively include, for example, the effect of cohesion (see Sec. 2.2 in
Sugiura et al. 2021). They demonstrated that, for $\phi_{\rm fri}\geq 70$ deg
with a spin-up timescale of $\lesssim$ a few days (defined as the elapsed time
of spin-up from the rotation period of 3.5 h to 3.0 h), an axisymmetric set of
surface landslides would occur and top-shaped rubble-pile bodies would be
correspondingly formed (panels (a)-(b) in Fig. 1).

The story of the landslide hypothesis, however, should not have ended here,
although Sugiura et al. (2021) did not study the fate of the surface materials
of the landslide. In this study, we first show that the surface materials are
distributed around the newly formed top-shaped central body, forming a
transient particle disk (panel (c) in Fig. 1). Here, the disk mass is
typically found to be about $\sim 10$wt.% of the top-shaped primary body.

In this study, following the numerical approach used in Sugiura et al. (2021),
we investigate a full long-term dynamical evolution – from the landslide to
the end of the disk evolution – to understand the fate of the surface
materials that are distributed around the newly formed top-shaped body (Fig.
1). We show that moons and an axisymmetric equatorial ridge are naturally and
inevitably formed around the top-shaped primary (panel (d) in Fig. 1).

In Sect. 2, we describe our numerical models. In Sect. 3, we show our
numerical results. Section 4 discusses the long-term orbital evolution between
a moon and the central primary due to the binary YORP (BYORP) effect, which
may lead to diverse systems of the top-shaped asteroids with and without a
companion moon(s) (panel (e) and panel (f) in Fig. 1). Finally, Section 5
summarizes this paper.

##  2 Numerical method

Our SPH simulations are the same as those in Sugiura et al. (2021), and we
refer the readers to their Sec. 2 for more details. Our simulations solve
hydrodynamic equations (the equations of continuity and motion) with self-
gravity and yield processes with the aid of a friction model that determines
the shear strength of granular material (Jutzi, 2015).

In a similar context, a Discrete Element Methods (DEMs) with cohesionless
particles were used in Walsh et al. (2008). Shape deformation using this
approach could be numerically affected due to the inconsistency between the
realistic particle size and that used in their simulations. In the SPH
approach, in contrast, a rubble-pile body was constructed by a continuum of
granular material and we explicitly set a specific value of the friction angle
of the material. We expect that the SPH method has the advantage of precise
control of bulk friction of rubble piles (see also Jutzi & Asphaug, 2015). We
also note that there is a difference in the process of the moon formation
between their study and ours; in Walsh et al. (2008), the moon formation
occurred via a step-by-step accumulation of discretely ejected particles,
while, in our study, the moon formation occurs because of particle disk
evolution.

The rubble-pile body initially had a uniform and spherical body with a radius
of $R_{\rm ini}=500$ m, which was then numerically spun up with the angular
acceleration of $8.954\times 10^{-10}$ rad s-2. We stopped the acceleration
when 1wt.% of the primary is ejected; this prevents further artificial
deformation after the landslide, and the arbitrary choice of 1wt.% does not
affect our results as long as the landslide occurs with a shorter timescale
(here $\sim 3$ h) than the spin-up timescale (here $\sim 30$ h). The density
of particles at uncompressed states is $\rho_{\rm ini}=1.19$ g cm-3. Particle
mass is given as $m=M_{\rm tot}/N_{\rm tot}$, where $M_{\rm
tot}\equiv(4/3)\pi\rho_{\rm ini}R_{\rm ini}^{3}$ and $N_{\rm tot}=25470$ is
the total number of the SPH particles.

In this study, following Sugiura et al. (2021), we specifically focused on a
set of parameters (i.e., $\phi_{\rm fri}=80$ deg; see panels (k)-(o) in their
Fig. 1), which successfully resulted in the formation of an axisymmetric top-
shaped figure as seen at asteroids Ryugu, Bennu, and Didymos. Lowering
$\phi_{\rm fri}$, for example, resulted in an elongated lemon-shaped primary
via internal deformation (see panels (a)-(e) and (f)-(j) in Fig. 1 of Sugiura
et al. 2021), which eventually led to distribute particles as disks with a
variety of initial disk mass. The dependences on the disk parameters (e.g.,
initial disk mass) and tidal parameters are studied in detail using an
independent numerical approach (1D fluid simulations) in our companion paper
(Madeira et al. submitted). Here, we used the direct and continuous approach –
from the landslide to the end of the disk evolution – and focused on studying
the details regarding the relevant consequences of the top-shape formation.

Below, we show that the surface materials of the primary are distributed in a
disk-like structure and such a particulate disk eventually forms moon(s)
around the top-shaped primary. We numerically distinguished member particles
of the primary as follows. We start from a randomly chosen particle and then
iteratively detect any particles within a critical distance in a bottom-up
fashion (Hyodo & Genda, 2018). This procedure continues until no more new
particles are detected. As a critical distance, we used 1.5 times the
smoothing length of an SPH particle to be consistent with Sugiura et al.
(2021), although this specific choice does not affect our results.

In this study, we continued our simulation up to $t\sim 25\times 10^{5}$ s
($\sim 700$ hours). This is much after the landslide had occurred; Sugiura et
al. (2021) stopped their simulations at $t\sim 1\times 10^{5}$ s ($\sim 30$
hours) when the top-shaped figure was formed. We may need to be careful about
the angular momentum (AM) conservation in SPH simulations, particularly for a
rotating system (such as disk evolution). In our calculations, after an
artificial spin-up of the primary (i.e., artificial AM increase in the system)
was stopped, the cumulative AM error was $\sim 2.6$% during and after the disk
formation until the end of the simulation (i.e., from $t\sim 1.5\times 10^{5}$
s to $t\sim 25\times 10^{5}$ s).

Although our SPH simulations contain some error in AM conservation,
independent numerical simulations of particle disk evolutions of the Earth’s
Moon formation (e.g., $N$-body simulations; e.g., Kokubo et al. 2000 and 1D
semi-analytical fluid simulations; e.g., Salmon & Canup 2012) used similar
initial disk conditions (e.g., in terms of disk to primary mass ratio) and
their results (e.g., resultant moon to primary mass ratio) were consistent
with those obtained here. This is because the gravitational disk evolution can
be well characterized by the mass ratio between the disk and the primary
(Kokubo et al., 2000; Salmon & Canup, 2012). Furthermore, using $N$-body
simulations, Hyodo et al. (2015) showed a scaling law for the mass of the
largest moon formed via disk spreading as a function of the initial disk mass
(see their Fig. 12) and it is generally consistent with that obtained in this
study. We note that, although these studies generally assumed particulate
disks, a more realistic Moon-forming disk would be vapor-rich and thus the
Earth’s Moon formation process would be more complex (e.g., Thompson &
Stevenson, 1988; Nakajima & Stevenson, 2014; Lock et al., 2018)

Figure 2: Overall disk evolution after the landslide. Gray particles represent
member particles of the top-shaped primary. Dark-gray particles represent disk
particles and member particles of gravitational clumps (including the formed
moon at later epochs). The black circle indicates the Roche limit of the
primary ($r_{\rm R}\sim 2.5R_{\rm c}$ where $R_{\rm c}=500$m is used).  Figure
3: Same as Fig. 2, but edge on views.

##  3 Numerical results: moon formation

Figure 4: Evolutions of the disk mass (left panel; not including the mass of
the largest moons) and the masses of the largest three moons (right panel). At
around $t\sim 1\times 10^{5}$s with a time-interval of $\Delta t\sim 10^{4}$s,
the landslide occurs and a particle disk is formed (left panel). The disk
spreads and starts to form gravitational clumps with masses larger than $\sim
0.001M_{\rm c}$ at $t\sim 5.5\times 10^{5}$s, where $M_{\rm c}$ is the mass of
the central primary. Two large moons are formed by the end of our simulations
($t=25\times 10^{5}$s; $\sim 700$h).

Figures 2 and 3 show an overall evolution of a particle disk including the
landslide. The landslide (i.e., disk formation) occurs with a very short
timescale ($\sim 10^{4}$ s $\simeq 3$ hours) compared to the disk evolution
timescale ($\sim 10^{6}$ $\simeq 300$ hours; see Fig. 4). Because of the
landslide, a top-shaped figure of the primary is formed (see more details in
Sugiura et al. 2021 and their Fig.1). The surface materials of the primary are
distributed around the primary, forming a particle disk. The initial mass of
the disk is $M_{\rm disk}\sim 20$wt.% of the top-shaped primary (Fig. 4). Most
of the disk mass is initially within the Roche limit of the primary, defined
as $r_{\rm R}\equiv 2.456(\rho_{\rm m}/\rho_{\rm c})^{-1/3}R_{\rm c}$ where
$\rho_{\rm m}$ and $\rho_{\rm c}$ are densities of the moon and the central
body. $R_{\rm c}$ is the radius of the central body.

The disk mass is large enough to become gravitationally unstable and the
spiral arm structures are formed as a result of gravitational instability.
This leads to an efficient angular momentum transfer and the outer disk
spreads further radially outward (Takeda & Ida, 2001), distributing the disk
materials beyond the Roche limit. Because of the angular momentum
conservation, the inner part of the disk spreads inward, resulting in re-
accretion onto the primary.

The materials scattered beyond the Roche limit start to coagulate via their
own gravity, forming gravitational clumps (Figs. 2 and 3). With time,
accretion among clumps proceeds, forming larger gravitational bodies (see also
the right panel of Fig. 4 for the mass evolution of the largest objects).
After $\sim 300$ hours since the landslide, large moons around the top-shaped
primary asteroid are formed. The total mass of the formed moons is $\sim
4$wt.% of the primary.

Here, two large moons are formed (Figs. 2, 3, and 4). The number of the large
moon is the result of the stochastic nature of the accretion processes that
can be seen in the 3D simulations, especially when the disk is massive to the
central body; the details of the stochastic nature were reported in
independent $N$-body simulations (Hyodo et al., 2015). This means that, if one
runs another simulation with the same parameters but with slightly different
initial positions of particles, the final outcome could be the formation of a
single moon. Thus, the number of large moons is not a decisive parameter that
characterizes the general outcome. More importantly, instead, regardless of
the number of moons, the total mass of the moons rarely changes and it
characterizes the outcome of specific choices of the disk parameters (e.g.,
initial disk mass; see Ida et al. 1997; Kokubo et al. 2000).

##  4 Discussion

###  4.1 Formation of equatorial ridge

Figure 5: Comparison of asteroid shapes between the initial figure ($t=0$h;
left panel) and after the disk evolution ($t\sim 700$h; right panel). A black
circle indicates a circle with a radius of 500 meters.

Our numerical simulations show another interesting feature as a result of disk
evolution around a top-shaped asteroid. Because the disk spreads both inward
and outward to conserve the angular momentum, the disk particles initially at
the inner region re-accrete onto the top-shaped primary. This has led to the
formation of an equatorial ridge.

Other potential physical mechanisms to form an equatorial ridge are discussed
in other literature. These are, for example, rotational reshaping and mass
movement (e.g., Walsh et al., 2012; Hirabayashi & Scheeres, 2015; Hirabayashi
et al., 2020), re-accumulation of debris after a catastrophic disruption of a
parent body (Michel et al., 2020), or selective accumulation of ejecta of
small impacts at rapidly rotating asteroid (Ikeya & Hirata, 2021).

Figure 5 shows an edge on views of the central rubble-pile body at $t=0$h
(left panel) and at the end of our simulation without disk particles and
formed moons ($t\sim 700$h; right panel). A significant shape change can be
seen from the initial spherical shape to a top-shaped figure with a prominent
axisymmetric equatorial ridge. This is a direct consequence of the re-
accretion of equatorial disk particles onto the top-shaped primary.

Many of the top-shaped asteroids and relatively spheroidal bodies in the near-
Earth region are reported to have an equatorial ridge (Benner et al., 2015).
JAXA’s Hayabusa2 mission revealed that asteroid Ryugu is a top-shaped asteroid
with an equatorial ridge (Watanabe et al., 2019). NASA’s OSIRIS-REx mission
also reported that asteroid Bennu has an equatorial ridge (Walsh et al.,
2019). Another top-shaped asteroid, Didymos, has also an equatorial ridge as
well as a moon, Dimorphos (Naidu et al., 2020).

We importantly note that the axis ratio of $c/a$ and the equatorial ridge seen
in Fig. 5 are smaller and much more prominent than the observations of, for
example, Ryugu, Bennu, and Didymos (Watanabe et al., 2019; Lauretta et al.,
2019; Naidu et al., 2020). Our SPH simulations would be limited in their
ability to precisely demonstrate the detailed accumulation processes of
particles in the equatorial ridge; because our SPH simulations assumed and
employed, for example, the same size and simple physical properties among all
particles. We also emphasize the importance of studying the long-term
geological evolution (e.g., including degradation) of the ridge. Continuous
micrometeoroid impacts and/or thermal fatigue would also change the ridge
shape. Therefore, further studies of the detailed ridge formation process as
well as the post-formation geological processes need to be done to fully
validate our ridge formation scenario via the landslides followed by the disk
evolution. We leave these points to later work.

###  4.2 On the diversity of top-shaped asteroids

Some top-shaped asteroids, e.g., Didymos, are found to have a companion
moon(s), while others do not (e.g., Ryugu and Bennu). In the discussion below,
we focus on asteroids that have (1) top-shaped figures and (2) equatorial
ridges, as background considerations. We, then, additionally consider the
existence and non-existence of a companion moon around the top-shaped primary.

Our results together with those of Sugiura et al. (2021) indicate that small
rubble-pile asteroids may inevitably experience a landslide due to spin-up by
the YORP effect or by other physical processes (e.g., re-accumulation, small
impact, and/or a close encounter with a planet). Depending on the acceleration
rate of spin-up as well as the effective friction angle, $\phi_{\rm fri}$, of
constituent particles of an asteroid, the resultant shape of asteroids via a
landslide would change (e.g., a lemon-shape or a top-shape; see Fig.1 of
Sugiura et al. 2021).

When the top-shaped primary was formed with the fast spin-up rate and
$\phi_{\rm fri}\gtrsim 70$ deg, Sugiura et al. (2021) numerically demonstrated
that about $\sim 10$wt.% of the primary is generally ejected. When the initial
disk mass is $\sim 10$wt.% or larger, the mass of the moon would be linearly
scaled with the disk mass (Ida et al., 1997; Kokubo et al., 2000), and thus a
change in the resultant disk mass within the same order of magnitude does not
significantly change the mass of the moon (would be on the order of $\sim
1$wt.% of the primary as observed for the top-shaped binary asteroids). In
contrast, when the effective friction angle is smaller than the above value,
the shape deformation became more gradual, and a lemon-shaped primary was
formed. In this case, the ejected mass tends to be smaller (e.g., $\sim 3$wt.%
and $\sim 5$wt.% for $\phi_{\rm fri}=40$ deg and $\phi_{\rm fri}=60$ deg,
respectively), although more detailed studies on throughout parameters may be
needed (see Fig. 1 of Sugiura et al., 2021). In these cases of smaller initial
disk masses, the mass of the resultant moon depends more strongly on the disk
mass and much smaller moons can be produced from the disk (Hyodo et al.,
2015).

Importantly, as a natural consequence of a landslide, a particle disk is
formed around the primary. Once a particle disk exists, the disk inevitably
spreads and forms moon(s) as long as the disk is massive enough. The final
mass and orbital configurations of the moon systems formed through the disk
spreading depend on the initial disk mass and tidal parameters. Such
dependencies are studied using 1D fluid simulations in our companion paper in
the context of the Didymos-Dimorphos system formation (Madeira et al.
submitted). In short, the direct consequence of the landslide would be the
system of a moon(s) around a central asteroid – regardless of the top-shaped
or the lemon-shaped central primary – with an equatorial ridge on the primary.

Then, the question now is can we remove the formed moon if the top-shaped
figures of Ryugu and/or Benuu are formed via landslides? This is because these
top-shaped asteroids today do not have a moon around them. Ejection via tidal
evolution would not be promising because its efficiency significantly
decreases as the distance between the primary and the moon becomes larger.

A potential dominant dynamical mechanism on the binary system may be the
binary YORP (BYORP) effect (e.g., Ćuk & Burns, 2005; Ćuk & Nesvorný, 2010;
McMahon & Scheeres, 2010; Jacobson & Scheeres, 2011). Although the timescale
of orbital separation to eject a companion moon via the BYORP effect can be as
small as $\sim 10^{5}$ years (see Sec.4.3 in Sugiura et al. 2021 and Ćuk &
Burns, 2005), its timescale and direction of the orbital separation (shrinking
or expanding) strongly depend on the details of the surface properties of the
asteroids and moons (Ćuk & Burns, 2005; McMahon & Scheeres, 2010).

We envision that such a complex dependence on the shapes and surface
properties can potentially lead to a diversity of the top-shaped asteroids
with and without a companion moon. The moons can be ejected by the BYORP
effect in some cases, while the dynamical evolution of the orbital separation
may not be efficient in other cases, remaining the system as binary asteroids
(Figure 1).

##  5 Summary

In this study, using the SPH simulations, we studied a continuous dynamical
sequence of events from the surface landslide that forms a top-shaped asteroid
and a particulate disk to disk evolution (Figure 1). We numerically
demonstrated that the particle disk formed by a surface landslide would
quickly spread and produce moons just outside the Roche limit of the top-
shaped primary (within $\sim 300$ hours). The mass of the moon is consistent
with what is observed around the top-shaped asteroids (on the order of $\sim
1$wt.% of the primary). We also demonstrated that an equatorial ridge would be
naturally formed because a fraction of the disk particles re-accrete
selectively on the equatorial region of the primary.

Tidal interaction as well as the binary YORP (BYORP) effect between moons and
the primary would change the orbital separation between them. The timescale
and direction of the orbital separation (shrinking or expanding) strongly
depend on the details of the surface properties of the moons and the primary.
This indicates that a long-term orbital evolution could produce diverse moon
systems around the top-shaped asteroids.

We envision that top-shaped Ryugu and Bennu could once have an ancient moon
that was later lost as a result of a successive moon’s orbital evolution.
Alternatively, other top-shaped asteroids today, such as Didymos, have a
companion moon. In these cases, no significant orbital evolution of the moon
may have occurred that would result in its loss.

Our study focused on the top-shaped asteroids. However, our results of the
consequences of a surface landslide – the moon formation and the equatorial
ridge formation – can be qualitatively (but not always quantitatively) applied
to any small rubble-pile asteroids near the breakup rotation period. Indeed,
Sugiura et al. (2021) showed that changing the effective friction angle and/or
the spin-up rate results in different patterns of the deformation mode with a
variety of disk formation (see their Fig. 1). Once a particle disk is formed,
the moon and the equatorial ridge would be naturally and inevitably formed.
The long-term evolution, then, would lead to a diverse asteroid system (see
also Ćuk, 2007; Jacobson & Scheeres, 2011; Ćuk et al., 2021).

Further development of the theoretical and modeling research, especially on
the BYORP effect, together with a better understanding of the surface and
particle properties of individual asteroids would be needed to further
constrain and validate the results presented in this study. Further studies on
the accumulation and post-formation geological processes of the equatorial
ridge are also demanded. Analysis of the return samples by JAXA’s Hayabusa2
and NASA’s OSIRIS-REx as well as data that would be obtained by NASA’s DART
and ESA’s HERA missions would help us to better understand the nature and
evolution of the top-shaped asteroids.

R.H. acknowledges the financial support of MEXT/JSPS KAKENHI (Grant Number
JP22K14091). R.H. also acknowledges JAXA’s International Top Young program.
K.S. acknowledges the financial support of JSPS, Japan KAKENHI Grant
(JP20K14536, JP20J01165). Numerical simulations in this work were carried out
on the Cray XC50 supercomputer at the Center for Computational Astrophysics,
National Astronomical Observatory of Japan.

## References

  * Becker et al. (2015) Becker, T. M., Howell, E. S., Nolan, M. C., et al. 2015, Icarus, 248, 499, doi: 10.1016/j.icarus.2014.10.048
  * Benner et al. (2015) Benner, L. A. M., Busch, M. W., Giorgini, J. D., Taylor, P. A., & Margot, J. L. 2015, in Asteroids IV, 165–182, doi: 10.2458/azu_uapress_9780816532131-ch009
  * Brozović et al. (2011) Brozović, M., Benner, L. A. M., Taylor, P. A., et al. 2011, Icarus, 216, 241, doi: 10.1016/j.icarus.2011.09.002
  * Ćuk (2007) Ćuk, M. 2007, ApJ, 659, L57, doi: 10.1086/516572
  * Ćuk & Burns (2005) Ćuk, M., & Burns, J. A. 2005, Icarus, 176, 418, doi: 10.1016/j.icarus.2005.02.001
  * Ćuk et al. (2021) Ćuk, M., Jacobson, S. A., & Walsh, K. J. 2021, , 2, 231, doi: 10.3847/PSJ/ac3093
  * Ćuk & Nesvorný (2010) Ćuk, M., & Nesvorný, D. 2010, Icarus, 207, 732, doi: 10.1016/j.icarus.2009.12.005
  * Harris et al. (2009) Harris, A. W., Fahnestock, E. G., & Pravec, P. 2009, Icarus, 199, 310, doi: 10.1016/j.icarus.2008.09.012
  * Hirabayashi & Scheeres (2015) Hirabayashi, M., & Scheeres, D. J. 2015, ApJ, 798, L8, doi: 10.1088/2041-8205/798/1/L8
  * Hirabayashi et al. (2020) Hirabayashi, M., Nakano, R., Tatsumi, E., et al. 2020, Icarus, 352, 113946, doi: 10.1016/j.icarus.2020.113946
  * Hyodo et al. (2016) Hyodo, R., Charnoz, S., Genda, H., & Ohtsuki, K. 2016, ApJ, 828, L8, doi: 10.3847/2041-8205/828/1/L8
  * Hyodo & Genda (2018) Hyodo, R., & Genda, H. 2018, ApJ, 856, L36, doi: 10.3847/2041-8213/aab7f0
  * Hyodo et al. (2015) Hyodo, R., Ohtsuki, K., & Takeda, T. 2015, ApJ, 799, 40, doi: 10.1088/0004-637X/799/1/40
  * Ida et al. (1997) Ida, S., Canup, R. M., & Stewart, G. R. 1997, Nature, 389, 353, doi: 10.1038/38669
  * Ikeya & Hirata (2021) Ikeya, R., & Hirata, N. 2021, Icarus, 367, 114590, doi: 10.1016/j.icarus.2021.114590
  * Jacobson & Scheeres (2011) Jacobson, S. A., & Scheeres, D. J. 2011, ApJ, 736, L19, doi: 10.1088/2041-8205/736/1/L19
  * Jutzi (2015) Jutzi, M. 2015, Planet. Space Sci., 107, 3, doi: 10.1016/j.pss.2014.09.012
  * Jutzi & Asphaug (2015) Jutzi, M., & Asphaug, E. 2015, Science, 348, 1355, doi: 10.1126/science.aaa4747
  * Kokubo et al. (2000) Kokubo, E., Ida, S., & Makino, J. 2000, Icarus, 148, 419, doi: 10.1006/icar.2000.6496
  * Lauretta et al. (2019) Lauretta, D. S., Dellagiustina, D. N., Bennett, C. A., et al. 2019, Nature, 568, 55, doi: 10.1038/s41586-019-1033-6
  * Lock et al. (2018) Lock, S. J., Stewart, S. T., Petaev, M. I., et al. 2018, Journal of Geophysical Research (Planets), 123, 910, doi: 10.1002/2017JE005333
  * Margot et al. (2015) Margot, J. L., Pravec, P., Taylor, P., Carry, B., & Jacobson, S. 2015, in Asteroids IV, 355–374, doi: 10.2458/azu_uapress_9780816532131-ch019
  * McMahon & Scheeres (2010) McMahon, J., & Scheeres, D. 2010, Celestial Mechanics and Dynamical Astronomy, 106, 261, doi: 10.1007/s10569-009-9247-9
  * Michel et al. (2020) Michel, P., Ballouz, R. L., Barnouin, O. S., et al. 2020, Nature Communications, 11, 2655, doi: 10.1038/s41467-020-16433-z
  * Naidu et al. (2020) Naidu, S. P., Benner, L. A. M., Brozovic, M., et al. 2020, Icarus, 348, 113777, doi: 10.1016/j.icarus.2020.113777
  * Nakajima & Stevenson (2014) Nakajima, M., & Stevenson, D. J. 2014, Icarus, 233, 259, doi: 10.1016/j.icarus.2014.01.008
  * Ostro et al. (2006) Ostro, S. J., Margot, J.-L., Benner, L. A. M., et al. 2006, Science, 314, 1276, doi: 10.1126/science.1133622
  * Rubincam (2000) Rubincam, D. P. 2000, Icarus, 148, 2, doi: 10.1006/icar.2000.6485
  * Salmon & Canup (2012) Salmon, J., & Canup, R. M. 2012, ApJ, 760, 83, doi: 10.1088/0004-637X/760/1/83
  * Sugiura et al. (2021) Sugiura, K., Kobayashi, H., Watanabe, S.-i., et al. 2021, Icarus, 365, 114505, doi: 10.1016/j.icarus.2021.114505
  * Takeda & Ida (2001) Takeda, T., & Ida, S. 2001, ApJ, 560, 514, doi: 10.1086/322406
  * Takeda & Ohtsuki (2009) Takeda, T., & Ohtsuki, K. 2009, Icarus, 202, 514, doi: 10.1016/j.icarus.2009.03.001
  * Thompson & Stevenson (1988) Thompson, C., & Stevenson, D. J. 1988, ApJ, 333, 452, doi: 10.1086/166760
  * Walsh & Jacobson (2015) Walsh, K. J., & Jacobson, S. A. 2015, in Asteroids IV, 375–393, doi: 10.2458/azu_uapress_9780816532131-ch020
  * Walsh et al. (2008) Walsh, K. J., Richardson, D. C., & Michel, P. 2008, Nature, 454, 188, doi: 10.1038/nature07078
  * Walsh et al. (2012) —. 2012, Icarus, 220, 514, doi: 10.1016/j.icarus.2012.04.029
  * Walsh et al. (2019) Walsh, K. J., Jawin, E. R., Ballouz, R. L., et al. 2019, Nature Geoscience, 12, 242, doi: 10.1038/s41561-019-0326-6
  * Watanabe et al. (2019) Watanabe, S., Hirabayashi, M., Hirata, N., et al. 2019, Science, 364, 268, doi: 10.1126/science.aav8032