# Navigating Fairness in AI Radiology: Concepts, Consequences, and Crucial
Considerations

Vasantha Kumar Venugopal Corresponding author<EMAIL_ADDRESS>CARPL.ai,
New Delhi, India  Abhishek Gupta Corresponding author:
<EMAIL_ADDRESS>CARPL.ai, New Delhi, India  Rohit Takhar
CARPL.ai, New Delhi, India  Charlene Liew Jin Yee  Duke-NUS Medical School,
Singapore  Catherine Jones  I-MED Radiology Network, Australia  Gilberto Szarf
Hospital Israelita Albert Einstein, Brazil

###### Abstract

Artificial Intelligence (AI) has significantly revolutionized radiology,
promising improved patient outcomes and streamlined processes. However, it’s
critical to ensure the fairness of AI models to prevent stealthy bias and
disparities from leading to unequal outcomes. This review discusses the
concept of fairness in AI, focusing on bias auditing using the Aequitas
toolkit, and its real-world implications in radiology, particularly in disease
screening scenarios.

Aequitas, an open-source bias audit toolkit, scrutinizes AI models’ decisions,
identifying hidden biases that may result in disparities across different
demographic groups and imaging equipment brands. This toolkit operates on
statistical theories, analyzing a large dataset to reveal a model’s fairness.
It excels in its versatility to handle various variables simultaneously,
especially in a field as diverse as radiology. The review explicates essential
fairness metrics: Equal and Proportional Parity, False Positive Rate Parity,
False Discovery Rate Parity, False Negative Rate Parity, and False Omission
Rate Parity. Each metric serves unique purposes and offers different insights.
We present hypothetical scenarios to demonstrate their relevance in disease
screening settings, and how disparities can lead to significant real-world
impacts.

##  1 Key Points

  1. 1.

The Aequitas toolkit is a powerful resource for auditing and mitigating bias,
providing a wide range of metrics to help dissect and understand the potential
biases lurking in AI models.

  2. 2.

Different fairness metrics serve distinct purposes in AI systems. In a disease
screening scenario, metrics such as Equal and Proportional Parity, False
Positive Rate Parity, False Discovery Rate Parity, False Negative Rate Parity,
and False Omission Rate Parity play significant roles

  3. 3.

Understanding each metric’s nuances is essential for ensuring a comprehensive
approach to bias mitigation.

##  2 Introduction

Artificial intelligence (AI) has staked its claim as a game-changer in
healthcare, with an extraordinary potential to transform various aspects,
including disease diagnosis, treatment, and management [AKB+19, Top19,
KLvG+17]. Radiology, a field inherently dependent on image interpretation, has
emerged as a prime candidate for the AI revolution [CKM+18, SWS17]. Indeed, AI
algorithms’ ability to process and learn from vast volumes of imaging data
offers an enticing prospect for enhancing radiologists’ diagnostic accuracy
and efficiency [UYS+21]. However, as we integrate AI deeper into radiology,
it’s essential to ensure these systems operate fairly, delivering consistent
performance across diverse patient groups. Unfortunately, this isn’t a given.

AI models can unwittingly perpetuate or even amplify existing biases, leading
to disparities in healthcare outcomes [OPVM19, RHH+18]. Such biases can
originate from various sources, including the data used to train these models
and the methods employed in their development [9]. It’s a pressing concern
that has garnered significant attention, prompting researchers to develop
methods and tools to detect and mitigate AI biases [RHH+18, PMA19]. This
review will delve into these critical areas, focusing on the concept of bias,
disparity, and fairness, the Aequitas toolkit - a key player in assessing AI
fairness, and the real-world impacts of these issues in radiology AI.

##  3 Disparity, Bias, and Fairness – Definitions and Explanations

Before diving into the mechanisms of bias in AI and its impact on radiology,
it is essential to establish a common understanding of some critical terms:
disparity, bias, and fairness. These terms, although sometimes used
interchangeably, hold distinct meanings in the context of AI and healthcare.

Disparity refers to differences, variations, or inequalities that occur
naturally or systematically within a population. In healthcare, these
disparities can emerge across various dimensions, including age, gender,
ethnicity, socioeconomic status, and more. Importantly, not all disparities
are unjust or undesirable. For instance, the higher incidence of breast cancer
in women compared to men is a biological disparity [MPMC19, SKH+18]. However,
if certain ethnic groups receive less accurate cancer screenings due to
algorithmic bias, this systematic disparity becomes a concern.

Bias, in the context of AI, refers to the systematic error introduced into the
model that leads to unfair treatment or outcomes for certain groups. This
could be due to several factors, such as skewed data representation, biased
labels, or the algorithm’s inherent limitations [Nel19]. This bias can be
unintentional and subtle, yet it can lead to unfair or unequal outcomes
[SKZM+21].

Fairness, a broad and multifaceted concept, generally refers to the state
where no group or individual is unduly disadvantaged or favored by an AI
system [CR18]. In healthcare AI, fairness can be understood as the equal and
impartial treatment of all individuals, regardless of their characteristics or
group affiliations [DHP+12]. Achieving fairness is challenging, especially
given the complex and multifaceted nature of healthcare data [ZKV+22].
However, it’s a pursuit worth striving for, considering the potential impacts
of unfair AI systems on people’s lives and health [RHH+18, PMA19].

Understanding these concepts forms the foundation for exploring bias in AI,
its detection, and its mitigation. In the following sections, we will delve
into these subjects more deeply, focusing on the Aequitas toolkit and its
relevance to radiology AI.

##  4 Aequitas: A Bias and Fairness Audit Toolkit

###  4.1 Aequitas Toolkit - What is it?

The Aequitas toolkit is an open-source bias audit tool, specifically created
to bring transparency and fairness into the machine learning models used in
various domains – healthcare, finance, criminal justice, and more. Developed
by the Center for Data Science and Public Policy at the University of Chicago,
Aequitas is a practical tool that invites us to take a closer look at the
fairness of our AI systems[SKH+18]. The name Aequitas, derived from Latin,
implies fairness and equality - and that’s precisely what this toolkit is
designed to promote [SKH+18].

Aequitas is developed with the aim of making machine learning models more
transparent, fair, and accountable [SKH+18]. It works by examining performance
data from a machine learning model and generating an array of metrics that
give insight into the model’s fairness across various subgroups.

###  4.2 Why Aequitas? What makes it the ideal tool for auditing bias in
radiology AI?

Firstly, it offers a comprehensive suite of fairness metrics, enabling a deep-
dive analysis of the different facets of fairness in a model. The broad range
of metrics allows for nuanced assessments of fairness, something that is
crucial in healthcare [OFG+20].

Secondly, Aequitas is designed to handle large and complex datasets, making it
perfect for the vast volumes of data involved in radiology AI. It assumes the
law of large numbers, meaning it works best when there’s a significant sample
size for each group you’re comparing.

Lastly, Aequitas doesn’t make assumptions about disease prevalence or
incidence among different groups, an essential feature when considering
fairness in disease screening scenarios. Instead, it looks at the model’s
predictions and outcomes across groups, allowing for a fair and unbiased
comparison.

##  5 Assessing Disparity and Parity

While disparity refers to the differences or imbalances in the distribution of
outcomes or predictions across different groups, its counterpart, parity
refers to the fairness or equity in those distributions.

###  5.1 Disparity Quantification

The aequitas bias toolkit employs disparate impact metrics for disparity
calculation. As a cardinal metric in disparity computation, particularly in
fairness analysis, disparate impact scrutinizes relative outcome or prediction
distributions across heterogeneous groups and detects possible biases. The
disparity is gauged relative to a reference group, given by:

$DisparityMeasure=\frac{OutcomeRateMetricsforspecificgroup}{OutcomeRateMetricsforreferencegroup}$

The Outcome Rate Metrics could be any precision metrics such as False Negative
Rate (FNR), False Omission Rate (FOR), Positive Predictive Value (PPV), etc. A
Disparity Measure of 1 indicates no disparity, whereas values diverging from 1
suggest the existence of disparity.

###  5.2 Reference Group Selection Strategy

The reference group selection is a vital aspect of the aequitas bias toolkit
analysis. It offers a standard against which other groups are juxtaposed to
detect disparities and assess fairness. Several strategies can be deployed to
select a reference group:

  1. 1.

Predominant Group as Reference: The most widespread group within the dataset
can serve as the reference group. It helps identify disparities or biases
relative to the dominant group.

  2. 2.

Demographically Balanced Reference: An averaged or balanced group,
representing the demographic or protected attribute under analysis, can serve
as the reference group to reflect the overall demographic constitution of the
study population

  3. 3.

External Benchmark: In specific instances, an external benchmark from existing
research or datasets can be the reference group.

  4. 4.

Customized Reference: Aequitas allows for reference group customization based
on specific criteria or research objectives.

The choice of the reference group should be aligned with the unique objectives
and considerations of the medical imaging analysis.

###  5.3 Parity

Parity denotes fairness or equity across different demographic or protected
groups. Parity achievement requires disparity minimization in the previously
mentioned metrics. Aequitas employs a threshold, the disparity intolerance,
representing the acceptable disparity level. When Disparity Measure values
surpass this threshold, it indicates unacceptable disparity, implying the need
for further inquiry and potential bias rectification. If the Disparity Measure
lies within the disparity intolerance range, the group is deemed in parity
with the reference group.

$\tau\leq\text{Disparity Measure}\leq\frac{1}{\tau}$ (1)

In equation 1, $\tau$ represents the disparity intolerance, within the range
[0,1]. By setting a specific limit, it offers a benchmark for evaluating
disparity magnitude, aiding in distinguishing between tolerable disparities
and those demanding intervention.

##  6 Fairness Metrics – Their Relevance in Disease Screening Scenario

Fairness metrics provide critical insights into the performance of AI models,
particularly in the context of disease screening. They act as evaluative
standards, ensuring that AI models are equitable across various patient
demographics [HCMD18].

###  6.1 Equal and Proportional Parity

These metrics ensure equal representation across diverse groups within a
dataset. These are crucial when you want to ensure representation among the
groups in your dataset. These metrics ensure that all groups have equal
chances of being flagged by the AI [Rie20, SKH+18]. However, in a disease
screening scenario, accurate disease identification is more important than
equal representation. So, while these metrics are interesting, they may not be
the main focus in a screening context.

Imagine we have a disease that is more prevalent in the older population. In
this case, if the AI system flags older individuals more often for further
testing, that is not necessarily a bias, but rather the algorithm correctly
identifying higher-risk individuals. However, if the AI system
disproportionately flags older individuals with a higher false positive rate,
this could be seen as a bias against the older population, leading to
unnecessary stress and potential harm from additional testing. Therefore,
while it’s important to consider equal representation, in a disease screening
scenario, it is even more crucial to ensure accurate disease identification
and avoid unnecessary harm.

###  6.2 False Positive Rate Parity

This metric evaluates the proportion of healthy individuals incorrectly
identified as having the disease within each demographic group [Raz23, 26-]. A
higher false positive rate in a particular group can lead to unnecessary
additional testing, inducing potential harm and increased healthcare costs.
FPR Parity ensures this proportion is the same across all groups [SKH+18].

This metric stipulates that each demographic group should experience the same
rate of false-positive errors. For instance, if we consider ethnicity as a
variable, false positive rate parity means that each ethnicity group—be it
Asian, African American, or Caucasian—should encounter an identical rate of
false positive errors.

The significance of this metric comes into focus particularly when the
consequences of false positives are severe and could lead to unnecessary
interventions or adverse outcomes. Hence, ensuring false positive rate parity
is crucial as it safeguards against disproportionately high false positive
rates in any single group.

In the context of tuberculosis (TB) screening using an AI-powered radiology
tool for immigration checks, the principle of False Positive Rate Parity is
quite significant. This concept demands equal false positive rates across all
demographic groups - for instance, various nationalities or ethnic groups
undergoing the immigration check.

Let’s consider a scenario where the AI system is being used to screen chest
X-rays of immigrants from various countries for signs of TB. A false positive
in this context would mean that an individual who does not have TB is
incorrectly identified by the AI system as having the disease. The consequence
of this error could be a delay in their immigration process, or even a denial
of immigration altogether. This could be seen as a punitive outcome, as it
unjustly impacts the individual’s prospects.

In this context, False Positive Rate Parity would mean that immigrants from
all countries have an equal likelihood of experiencing a false positive. This
ensures that no nationality or ethnic group is unduly affected by the errors
of the AI system, maintaining fairness in the immigration process.

The False Positive Rate is calculated as:

$FPR=\frac{FP}{FP+TN}$

where:

FP stands for False Positives (the number of negative instances incorrectly
classified as positive) & TN stands for True Negatives (the number of negative
instances correctly classified as negative)

###  6.3 False Discovery Rate Parity

This metric quantifies the proportion of false positives within the flagged
cases for each group [ZVRG16, SKH+18]. A higher false discovery rate indicates
a larger share of the flagged cases in a group being false positives, leading
to unnecessary medical interventions and patient anxiety.

This criterion asserts that every demographic group should experience an
equivalent rate of false discovery errors. For example, when analyzing the
factor of ethnicity, false discovery rate parity means all ethnic groups
should have the same rate of false discovery errors.

The significance of this criterion similar to FPR comes to the fore especially
when the outcome of a false discovery could lead to punitive consequences for
individuals.

False Discovery Rate Parity is another critical principle in the same scenario
of TB screening for immigration checks. This principle demands an equal rate
of false discovery errors across all demographic groups. In this context, a
false discovery would mean that among the individuals flagged by the AI system
as having TB (i.e., the “discovered” cases), an equal proportion across all
nationalities are actually false positives.

For instance, let’s assume the AI system flags 100 immigrants from Country A
and 50 immigrants from Country B as potential TB patients. If 30 out of the
100 flagged cases from Country A are false positives, then to maintain False
Discovery Rate Parity, the number of false positives from the 100 flagged
cases from Country B should also be around 15. This ensures that no
nationality bears an undue burden of false discoveries, maintaining fairness
in the immigration process.

The false discovery rate is the ratio of the number of false positive results
to the number of total predicted positive test results.

$FDR=\frac{FP}{PP}=\frac{FP}{FP+TP}$

where: FP stands for False Positives (the number of negative instances
incorrectly classified as positive), PP stands for Predicted Positives (the
number of positive predictions) & TP stands for True Positives (the number of
positive instances correctly classified as positive)

###  6.4 False Negative Rate Parity

This metric is of particular importance in disease screening scenarios. A
higher false negative rate indicates that the model overlooks more actual
disease cases in one group than in others [Raz23, 26-]. This bias could lead
to delayed treatment and poorer health outcomes for these patients.

In the context of radiology AI, False Negative Rate Parity is a significant
metric that calls for equal false negative rates across all identified groups
[OFG+20]. Let’s consider the example of gender groups in the patient
population.

Visualize a scenario where an AI-driven radiology tool is deployed for lung
cancer screening. A false negative in this context indicates that a patient
who indeed has lung cancer is inaccurately identified by the AI system as
cancer-free. This error can have serious implications as it may delay
necessary treatment for the patient, possibly leading to detrimental health
outcomes.

Let’s take gender groups as an example - men, women, and non-binary
individuals. If False Negative Rate Parity is achieved, it indicates that all
these gender groups have an equal likelihood of encountering a false negative
error. This is especially important when the role of the screening tool is to
assist in early disease detection and subsequent provision of treatment.

In essence, if the AI system is intended to help identify patients who require
medical intervention, such as lung cancer treatment, it becomes paramount to
ensure that it’s not missing patients from certain gender groups
disproportionately. Adherence to the principle of False Negative Rate Parity
can ensure that the AI system is equally effective across all gender groups,
thereby promoting fairness and reducing health disparities.

$FNR=\frac{FN}{FN+TP}$

where:

FN stands for False Positives (the number of positive instances incorrectly
classified as negative) & TP stands for True Positives (the number of positive
instances correctly classified as positive)

###  6.5 False Omission Rate Parity

False Omission Rate - FOR is the fraction of false negatives of a group within
the predicted negative of the group. It represents the proportion of actual
positive cases within the predicted negative cases for each group [SKH+18]. A
higher FOR suggests that a larger percentage of the non-flagged cases in a
group are actual positive cases, indicating missed opportunities for early
intervention [OFG+20] [ZVRG16].

Let’s understand this with an example. Consider a radiology AI tool tasked
with screening for osteoporosis, a condition that can occur across all age
groups, although it’s more common in older adults. In this context, a false
omission error happens when the AI tool fails to flag an individual who
actually has early signs of osteoporosis.

Let’s say, for instance, a middle-aged individual, who is at the onset of
osteoporosis, goes for a routine health checkup. The AI system, however, fails
to flag the early signs, and the individual is told they don’t need any
further osteoporosis-related intervention. This is a false omission. Here, the
person in question is actually in need of further examination and treatment
but has been erroneously classified as not needing any.

Now, if we apply the concept of False Omission Rate Parity, it would mean that
the AI system should have an equal probability of making this mistake across
all age groups – young adults (18-35), middle-aged adults (36-60), and older
adults (61 and above). This is crucial to prevent the AI tool from
disproportionately missing disease cases in certain age groups, thus ensuring
fairness in disease detection across all age groups.

By striving for False Omission Rate Parity, we can promote an equitable
allocation of healthcare resources and can better ensure that early signs of
diseases are not overlooked disproportionately in any age group.

$FOR=\frac{FN}{PN}=\frac{FN}{FN+TN}$

where: FN stands for False Negatives (the number of positive instances
incorrectly classified as negative), PN stands for Predicted Negative (the
number of Negative predictions) & TN stands for True Negatives (the number of
negative instances correctly classified as negative)

##  7 Hypothetical Scenarios to Demonstrate the Impact of Bias in Screening
Settings

Let’s imagine two hypothetical scenarios in medical screening settings where
we’ll observe the impact of bias. Here, we’ll bring the concepts of False
Positive Rate (FPR), False Discovery Rate (FDR), False Negative Rate (FNR),
and False Omission Rate (FOR) to life.

###  7.1 Scenario 1 to demonstrate the impact of FPR disparity: Tuberculosis
Screening for Australian Visa Applications

Australia’s academic appeal attracts a diverse range of international students
every year. In the 2021-22 period, the country received half a million student
visa applications. The applicants came from all over the globe, but the
majority hailed from five nations: China, India, the United Kingdom, the
United States, and Vietnam.

Figure 1: Hypothetical scenario to demonstrate the impact of the FPR disparity
using a TB Immigration screening scenario. The assumed FPR disparity leads to
a disproportionately huge penalty for the Indian immigrants

With 130,000 applications, Chinese students led the pack, making up 26% of the
total. Indian students followed closely with 110,000 applications, accounting
for 22% of all applications. The UK, the US, and Vietnam also made significant
contributions with 60,000 (12%), 50,000 (10%), and 40,000 (8%) applications,
respectively [Wor22] as shown in Figure 1.

Now, let’s introduce our AI-powered tuberculosis screening tool into this
immigration process. Let’s also assume that the tool has a sensitivity of 97%
and a specificity of 96% with a False Positive Rate (FPR) disparity towards
Indian applicants, with a disparity ratio of 1.5. This bias can significantly
affect the number of false positives among Indian applicants.

To illustrate this, let’s consider the TB prevalence in each of these
countries. China’s prevalence is 100 per 100,000 people, India’s is double
that at 200, the UK and the US have the lowest prevalence at 10 and 9
respectively, and Vietnam’s prevalence mirrors China’s at 100.

In the majority of cases, the tool performs well. For instance, in China, out
of the 130 expected cases of TB, the tool correctly identifies 126 but misses
4. It also falsely identifies 5,200 healthy individuals as having TB.
Similarly, in the UK and the US, the false positives are limited to 2,400 and
2,000, respectively.

However, the picture changes when we turn to India. Here, the expected cases
of TB rise to 220 in line with the high prevalence. The AI identifies 213
correctly but misses 7. The bias, though, results in a staggering 16,500
healthy individuals being flagged as having TB. That’s more than triple the
false positives in China, despite India having fewer applicants and only
double the TB prevalence.

This scenario highlights the substantial impact of bias in AI-powered disease
screening tools. The bias against Indian applicants can lead to a significant
increase in false positives, causing unnecessary stress and leading to a
higher burden of additional testing. Such disparities underscore the
importance of fairness in AI systems to prevent such disproportionate impacts.

###  7.2 Scenario 2 to demonstrate the impact of FNR disparity: Lung Cancer
Screening in Singapore

Singapore also has a diverse population, with Malays comprising about 15% of
the residents. With the introduction of an AI tool for lung cancer screening,
the hope was to ensure early detection and treatment across all ethnic groups.
However, let’s assume that our AI tool exhibits a False Negative disparity
with a 1.6 disparity ratio towards the Malay population.

With 100,000 individuals undergoing screening, about 2,000 people are expected
to have lung cancer given the 2% prevalence rate observed in global studies.

Figure 2: Flow chart describing the hypothetical FNR disparity scenario among
the Singaporean population

Let’s first consider the Chinese population, which represents 75.9% of the
total, or 75,900 individuals. With a disease prevalence of 2%, about 1,518
Chinese individuals would have lung cancer as shown in Figure 2. Given the AI
tool’s sensitivity of 98%, it would correctly identify about 1,487 of these
cases, missing around 31 cases. [18-22]

Now, let’s consider the Malays, who represent 15% of the total or 15,000
individuals. With a disease prevalence of 2%, about 300 Malay individuals
would have lung cancer. However, due to the AI tool’s bias, the effective
sensitivity for Malays drops to 61.25% (98%/1.6). This would mean the tool
correctly identifies only about 184 of these cases, missing around 116 cases
as in Figure 3.

Figure 3: Hypothetical scenario to demonstrate the impact of the FNR disparity
using a Lung cancer screening scenario. The assumed FNR disparity leads to a
disproportionately high ratio of missed cancers among Malay patients

Comparing the numbers, you can see the stark difference. Despite the Malay
population being about a fifth of the Chinese population, the number of missed
cases is nearly four times higher due to the bias in the AI tool. The impact
is significant, demonstrating how bias can lead to disproportionate negative
outcomes for certain groups.

##  8 Mitigation of bias

While a detailed discussion of bias mitigation strategies is beyond the scope
of this article, we want to present a broad overview of some recommended
strategies. A multi-pronged strategy is required for mitigating bias in AI.
Representative samples in training datasets must be ensured, AI learning
algorithms should be adjusted for bias reduction, and post-training
modifications to the model’s decisions need to be made. Transparency and
explainability in the AI decision-making process facilitate the pinpointing
and correcting of biases. Regular audits, such as those enabled by Aequitas,
serve as a vital tool for monitoring fairness over time.

##  9 Conclusion

Each of these metrics defined in the Aequitas package offers a unique
perspective on fairness, highlighting potential biases in AI models[SKH+18,
BDH+18]. In the context of disease screening, these metrics ensure that the
models are not only accurate on average but are also fair across all groups
[FKZ+22]. The hypothetical scenarios presented herein elucidate the profound
effects of compounded bias, underlining the importance of considering
intersectionality in our audits. Our examination demonstrated that bias, left
unchecked, could significantly skew diagnostic outcomes, potentially
exacerbating existing health disparities.

Nevertheless, it is important to remember that fairness is context-
specific[36]. The interpretation of these metrics must consider the specific
characteristics of the disease, the population, and the overall screening
strategy [MLCW23]. Therefore, the use of these metrics is a step towards a
more equitable healthcare system, where all patients, irrespective of their
demographics, receive accurate and timely disease detection [RKF+22, Kah22,
OFG+20].

## References

  * [18-22] Study visa statistics, 2022.
  * [26-] fairness: Fairness, accountability, and transparency in machine learning.
  * [AKB+19] Diego Ardila, Atilla P. Kiraly, Sujeeth Bharadwaj, Bokyung Choi, Joshua J. Reicher, Lily Peng, Daniel Tse, Mozziyar Etemadi, Wenxing Ye, Greg Corrado, David P. Naidich, and Shravya Shetty.  End-to-end lung cancer screening with three-dimensional deep learning on low-dose chest computed tomography.  Nature Medicine, 25(6):954–961, May 2019.
  * [BDH+18] Rachel K. E. Bellamy, Kuntal Dey, Michael Hind, Samuel C. Hoffman, Stephanie Houde, Kalapriya Kannan, Pranay Lohia, Jacquelyn Martino, Sameep Mehta, Aleksandra Mojsilovic, Seema Nagar, Karthikeyan Natesan Ramamurthy, John Richards, Diptikalyan Saha, Prasanna Sattigeri, Moninder Singh, Kush R. Varshney, and Yunfeng Zhang.  Ai fairness 360: An extensible toolkit for detecting, understanding, and mitigating unwanted algorithmic bias, 2018.
  * [CKM+18] Garry Choy, Omid Khalilzadeh, Mark Michalski, Synho Do, Anthony E. Samir, Oleg S. Pianykh, J. Raymond Geis, Pari V. Pandharipande, James A. Brink, and Keith J. Dreyer.  Current applications and future impact of machine learning in radiology.  Radiology, 288(2):318–328, August 2018.
  * [CR18] Alexandra Chouldechova and Aaron Roth.  The frontiers of fairness in machine learning, 2018.
  * [DHP+12] Cynthia Dwork, Moritz Hardt, Toniann Pitassi, Omer Reingold, and Richard Zemel.  Fairness through awareness.  In Proceedings of the 3rd Innovations in Theoretical Computer Science Conference. ACM, January 2012.
  * [FKZ+22] Shahriar Faghani, Bardia Khosravi, Kuan Zhang, Mana Moassefi, Jaidip Manikrao Jagtap, Fred Nugen, Sanaz Vahdati, Shiba P. Kuanar, Seyed Moein Rassoulinejad-Mousavi, Yashbir Singh, Diana V. Vera Garcia, Pouria Rouzrokh, and Bradley J. Erickson.  Mitigating bias in radiology machine learning: 3. performance metrics.  Radiology: Artificial Intelligence, 4(5), September 2022.
  * [HCMD18] J. Henry Hinnefeld, Peter Cooman, Nat Mammo, and Rupert Deese.  Evaluating fairness metrics in the presence of dataset bias, 2018.
  * [Kah22] Charles E. Kahn.  Hitting the mark: Reducing bias in AI systems.  Radiology: Artificial Intelligence, 4(5), September 2022.
  * [KLvG+17] Thijs Kooi, Geert Litjens, Bram van Ginneken, Albert Gubern-Mérida, Clara I. Sánchez, Ritse Mann, Ard den Heeten, and Nico Karssemeijer.  Large scale deep learning for computer aided detection of mammographic lesions.  Medical Image Analysis, 35:303–312, January 2017.
  * [MLCW23] Amarachi B. Mbakwe, Ismini Lourentzou, Leo Anthony Celi, and Joy T. Wu.  Fairness metrics for health AI: we have a long way to go.  eBioMedicine, 90:104525, April 2023.
  * [MPMC19] G. McCartney, F. Popham, R. McMaster, and A. Cumbers.  Defining health and health inequalities.  Public Health, 172:22–30, July 2019.
  * [Nel19] Gregory S Nelson.  Bias in artificial intelligence.  North Carolina medical journal, 80(4):220–222, 2019.
  * [OFG+20] Luis Oala, Jana Fehr, Luca Gilli, Pradeep Balachandran, Alixandro Werneck Leite, Saul Calderon-Ramirez, Danny Xie Li, Gabriel Nobis, Erick Alejandro Muñoz Alvarado, Giovanna Jaramillo-Gutierrez, Christian Matek, Arun Shroff, Ferath Kherif, Bruno Sanguinetti, and Thomas Wiegand.  Ml4h auditing: From paper to practice.  In Emily Alsentzer, Matthew B. A. McDermott, Fabian Falck, Suproteem K. Sarkar, Subhrajit Roy, and Stephanie L. Hyland, editors, Proceedings of the Machine Learning for Health NeurIPS Workshop, volume 136 of Proceedings of Machine Learning Research, pages 280–317. PMLR, 11 Dec 2020.
  * [OPVM19] Ziad Obermeyer, Brian Powers, Christine Vogeli, and Sendhil Mullainathan.  Dissecting racial bias in an algorithm used to manage the health of populations.  Science, 366(6464):447–453, October 2019.
  * [PMA19] Trishan Panch, Heather Mattie, and Rifat Atun.  Artificial intelligence and algorithmic bias: implications for health systems.  Journal of Global Health, 9(2), November 2019.
  * [Raz23] Shaina Raza.  Auditing icu readmission rates in an clinical database: An analysis of risk factors and clinical outcomes.  arXiv preprint arXiv:2304.05986, 2023.
  * [RHH+18] Alvin Rajkomar, Michaela Hardt, Michael D. Howell, Greg Corrado, and Marshall H. Chin.  Ensuring fairness in machine learning to advance health equity.  Annals of Internal Medicine, 169(12):866, December 2018.
  * [Rie20] Laura Rieger.  Explainability for neural networks.  2020\.
  * [RKF+22] Pouria Rouzrokh, Bardia Khosravi, Shahriar Faghani, Mana Moassefi, Diana V. Vera Garcia, Yashbir Singh, Kuan Zhang, Gian Marco Conte, and Bradley J. Erickson.  Mitigating bias in radiology machine learning: 1. data handling.  Radiology: Artificial Intelligence, 4(5), September 2022.
  * [SKH+18] Pedro Saleiro, Benedict Kuester, Loren Hinkson, Jesse London, Abby Stevens, Ari Anisfeld, Kit T. Rodolfa, and Rayid Ghani.  Aequitas: A bias and fairness audit toolkit, 2018.
  * [SKZM+21] Laleh Seyyed-Kalantari, Haoran Zhang, Matthew B. A. McDermott, Irene Y. Chen, and Marzyeh Ghassemi.  Underdiagnosis bias of artificial intelligence algorithms applied to chest radiographs in under-served patient populations.  Nature Medicine, 27(12):2176–2182, December 2021.
  * [SWS17] Dinggang Shen, Guorong Wu, and Heung-Il Suk.  Deep learning in medical image analysis.  Annual Review of Biomedical Engineering, 19(1):221–248, June 2017\.
  * [Top19] Eric J. Topol.  High-performance medicine: the convergence of human and artificial intelligence.  Nature Medicine, 25(1):44–56, January 2019.
  * [UYS+21] Daiju Ueda, Akira Yamamoto, Akitoshi Shimazaki, Shannon Leigh Walston, Toshimasa Matsumoto, Nobuhiro Izumi, Takuma Tsukioka, Hiroaki Komatsu, Hidetoshi Inoue, Daijiro Kabata, Noritoshi Nishiyama, and Yukio Miki.  Artificial intelligence-supported lung cancer detection by multi-institutional readers with multi-vendor chest radiographs: a retrospective clinical validation study.  BMC Cancer, 21(1), October 2021.
  * [Wor22] World Health Organization.  Global tuberculosis report, 2022.
  * [ZKV+22] Kuan Zhang, Bardia Khosravi, Sanaz Vahdati, Shahriar Faghani, Fred Nugen, Seyed Moein Rassoulinejad-Mousavi, Mana Moassefi, Jaidip Manikrao M. Jagtap, Yashbir Singh, Pouria Rouzrokh, and Bradley J. Erickson.  Mitigating bias in radiology machine learning: 2. model development.  Radiology: Artificial Intelligence, 4(5), September 2022.
  * [ZVRG16] Muhammad Bilal Zafar, Isabel Valera, Manuel Gomez Rodriguez, and Krishna P. Gummadi.  Fairness beyond disparate treatment and amp; disparate impact: Learning classification without disparate mistreatment.  2016\.


# Short-term memory by transient oscillatory dynamics in recurrent neural
networks

Kohei Ichikawa  Kunihiko Kaneko<EMAIL_ADDRESS>also, at
Research Center for Complex Systems Biology, Universal Biology Institute,
University of Tokyo.  Graduate School of Arts and Sciences, The University of
Tokyo, 3-8-1 Komaba, Meguro-ku, Tokyo 153-8902, Japan.

###### Abstract

Despite the significance of short-term memory in cognitive function, the
process of encoding and sustaining the input information in neural activity
dynamics remains elusive. Herein, we unveiled the significance of transient
neural dynamics to short-term memory. By training recurrent neural networks to
short-term memory tasks and analyzing the dynamics, the characteristics of the
short-term memory mechanism were obtained in which the input information was
encoded in the amplitude of transient oscillations, rather than the stationary
neural activities. This transient trajectory was attracted to a slow manifold,
which permitted the discarding of irrelevant information. Additionally, we
investigated the process by which the dynamics acquire robustness to noise. In
this transient oscillation, the robustness to noise was obtained by a strong
contraction of the neural states after perturbation onto the manifold. This
mechanism works for several neural network models and tasks, which implies its
relevance to neural information processing in general.

††preprint: APS/123-QED

##  I Introduction

Short-term memory is essential for our cognitive activitiesJonides  _et al._
(2008). Once an external signal is applied to an internal neuron, its
information is encoded therein and saved for a certain time period. Some form
of sustained neural activity is therefore necessaryColby  _et al._ (1996). A
conventional theory of working memory is based on the representation of each
memory item as a different attractor of the considered systemFunahashi  _et
al._ (1989); Fuster and Alexander (1971); Wang (2001). This theory is usually
referred to as “memories as attractors.” However, this is not practical when a
large number of objects has to be memorized, because it is difficult to
prepare so many attractors. Additionally, if there is a need to memorize
continuous information (e.g., length, frequency, amplitude of inputs), it
becomes more difficult to form a memory using this attractor. Furthermore,
extensive experimental reports suggest that while short-term memory is
maintained, neural activities are not constant; they instead continue to
change over timeHaider and McCormick (2009); Baeg  _et al._ (2003). Indeed,
the possibility of achieving short-term memory through sustained transient
dynamics has recently been discussedBondanelli and Ostojic (2020); Druckmann
and Chklovskii (2012); Goldman (2009); Orhan and Pitkow (2020); Orhan and Ma
(2019). Additionally, the general relevance of transient activities to neural
information processing has been discussedRabinovich  _et al._ (2008).

It is unclear, however, what form of transient neural activity can afford
short-term memory. In contrast to established studies on memories as
attractors, the way external inputs are encoded into transient dynamics is not
well explored. Furthermore, in contrast to stationary neural activities at
attractors, the way information is sustained in the transient process with
time-varying neural activities remains elusive. Memory must be robust under
noise to inputs or internal neural activities. From the view of memories as
attractors, the stability of memory is supported as the state is attracted to
the attractors (i.e., stationary states) after perturbationSussillo and Barak
(2013); Maheswaranathan  _et al._ (2019), whereas the robustness in the
transient dynamics is not well-explored. One must investigate how robustness
to noise is achieved to understand the transient dynamics that support short-
term memory.

To investigate the neural dynamics of short-term memory, we adopted a
recurrent neural network (RNN) trained to solve a task that requires
memorizing the input information for a given time span. We provided a task to
compare two interspaced signals input with some time interval. The RNN is
required to determine which of the two subsequent signals has the larger
continuous characteristics (e.g., frequencies of the periodic signals)Romo
_et al._ (1999). Here, to solve the comparison task, it is necessary for the
neural dynamics to maintain the information of the first signal until the
second signal input arrives. Thus, short-term memory is needed to solve this
comparison task. In this RNN, the neural dynamics consist of activities
entailing a large number of neurons connected via synapses whose weights are
adjusted by the standard back-propagation methodRumelhart  _et al._ (1986);
Werbos (1990) to solve this task.

After the network is trained to solve the task successfully, we analyze the
generated neural dynamics to uncover how the input information is encoded and
memorized according to the dynamics of neural activities. We find that the
neural activities exhibit a transient oscillation that endures for the time
span between the first and second signals. Subsequently, we analyze how this
oscillatory neural activity encodes the input information, provides short-term
memory, and solves the task. We uncover that the information of the first
signal (e.g., its frequency) is encoded and memorized by the amplitude of the
transient oscillation. Then, we determine if and how this memory by the
amplitude of transient oscillation is robust to noise using dynamical systems
theory.

The remainder of this paper is organized as follows. In the next section, we
introduce the task to compare the frequencies of two signals input with a
given time interval as well as the RNN model used to solve it. After
demonstrating that the trained RNN can solve the task, we analyze how the
memory of the first input is maintained. We demonstrate that neural activities
during the time span, while memory is maintained, exhibit a transient
oscillation. Here, the input-signal information required to solve the task is
encoded by the amplitude of this transient oscillation, whereas the irrelevant
information in the input signal that is unrelated to solving the task is
discarded. Then, the robustness of the memory to noise is analyzed. We also
confirm that the mechanism for short-term memory presented in this study is
valid for several different comparison tasks as well as in neural network
models that take biological features into the account, such as excitatory-
inhibitory balance or sparseness. Finally, the possible relevance of the
presented mechanism for short-term memory to biological neural dynamics is
discussed.

##  II model

###  II.1 Short-term memory task

As a task that requires short-term memory, we studied the frequency comparison
task illustrated in Fig.1, which is commonly adopted in the field of
neuroscienceRomo  _et al._ (1999). In this task, the input signals consist of
the first signal, a delay period, and the second signal. The objective of the
task is to determine which frequency is higher: that of the first or the
second signal.

Figure 1:  Schematic diagram of frequency comparison task. The input consists
of the first signal, a delay period during which there is no input signal, and
the second signal. Both the first and second signals are represented by noisy
sine waves. The task is to determine which of the frequencies of the first and
second signals is higher. The frequencies of two signals satisfy
$1\leq\omega_{1},\omega_{2}\leq 5$. In addition, only in the training phase
the difference of two signals satisfy $\ |\omega_{1}-\omega_{2}|\geq 1$.

Specifically, the first and second signals are chosen as noisy sine waves
following $u_{1,2}(t)=\sin(\omega_{1,2}t+\phi)+\eta_{1,2}(t)$, where
$\phi_{1,2}$ represents the phase of the signals, and $\eta_{1,2}(t)$ is a
random Gaussian variable with average zero and standard deviation 0.05. During
the delay period, there is no input signal. The duration of the first(second)
signal, $T_{s}(T_{s}^{\prime})$ and the delay period, $T_{d}$, vary with each
sample in the training phase. At $T_{s}$, the delay period starts, and at
$T_{f}=T_{s}+T_{d}$, the delay period ends. At
$T_{a}=T_{s}+T_{d}+T_{s}^{\prime}$, the second signal completes, and the
neural network is asked whether the first or second has the higher frequency.
Specifically, $T_{s}$ is homogeneously distributed as $T_{s}\in[13,17]$, and
$T_{d}$ is homogeneously distributed as $T_{d}\in[25,35]$ throughout the
samples in the training phase(see 1). During the test phase, to answer the
task, $T_{s}$ and $T_{d}$ are fixed at $15$ and $30$, respectively. These
conditions are fixed unless otherwise mentioned.

condition | training phase | test phase
---|---|---
Length of the input signal | $13\leq T_{s}\leq 17$ | $T_{s}=15$
Length of the delay period | $25\leq T_{d}\leq 35$ | $T_{d}=30$
Frequencies of signals | $1\leq\omega_{1},\omega_{2}\leq 5$ | $1\leq\omega_{1},\omega_{2}\leq 5$
Difference of $\omega_{1},\omega_{2}$ | $\ |\omega_{1}-\omega_{2}|\geq 1$ | $\ |\omega_{1}-\omega_{2}|\geq 0$
Table 1:  Conditions of the short-term memory task.

###  II.2 Recurrent Neural Network

In this study, a standard modelRivkind and Barak (2017) is adopted for neural
activity dynamics as expressed in the equation

$\dot{x}_{i}=-x_{i}+\sum_{j=1}^{N}J_{ij}\tanh(x_{j})+W_{i}^{\rm in}u,$ (1)

where $x_{i}$ represents the neural activity state of neuron $i$ (or $(1+\tanh
x_{i})/2$ can be correlated with the firing rate)111Because $x_{i}$ takes both
positive and negative values, it does not represent the neural activity per
se. Instead, in the context of neuroscience, $(1+\tanh(x_{i}))/2$ is often
regarded as the firing rate of the neuron.. Each neuron is recurrently
connected to the others, where $J_{ij}$ represents the strength of the
connection from neuron $j$ to $i$. Furthermore, $u_{i}$ represents the input
signal, which is defined by a cognitive task as described above and is
projected onto the neurons by the input connection, $W^{\rm in}_{i}$. The
number of neurons, $N$, is set at 256. The initial state of neuron
$x_{i}(t=0;i=1,...,N)$ is set to be a Gaussian random variable with average
zero and standard deviation of 0.1. In this study, the Euler method is applied
to simulate Eq. (1), and the discretized dynamics are subsequently calculated.
A discretized Eq.(1) is adopted here as
$x_{i}(t+1)=(1-\alpha)x_{i}(t)+\alpha(\sum_{j=1}^{N}J_{ij}\tanh(x_{j}(t))+W_{i}^{in}u(t))$.
The time width, $\alpha$, for discretization is set to 0.25; nonetheless, the
choice of this specific value is not essential to the results.

The output of the RNN is determined by the weighted sum of the neural states,
as stated in the equation,

$z_{i}(t)=\sum_{j=1}^{N}W_{ij}^{out}x_{j}.$ (2)

$\bf z$ is chosen to be a two-component vector, and $z_{1}-z_{2}$ represents
the logit of the probability that the first signal has a higher
frequencyBishop (2006). The calculation from ${\bf z}$ to the probability is
performed using the softmax function, ${\rm
Softmax}(z_{i})=e^{z_{i}}/(e^{z_{1}}+e^{z_{2}})$222The relationship between
logit $z_{1}-z_{2}$ and the probability calculated by the softmax function is
as follows. By definition of the softmax function, we have
$p_{1}=e^{z_{1}}/(e^{z_{1}}+e^{z_{2}})=1/(1+e^{-(z_{1}-z_{2})})$. Then, by the
definition of logit function, we get the logit of $p_{1}$ as
$\log(p_{1}/(1-p_{1}))=\log(e^{z_{1}-z_{2}})=z_{1}-z_{2}$. From the above
explanation, the output of the softmax function can be treated as a
probability because $p_{1}$ and $p_{2}$ satisfy $0\leq p_{1},p_{2}\leq 1$, and
$p_{1}+p_{2}=e^{z_{1}}/(e^{z_{1}}+e^{z_{2}})+e^{z_{2}}/(e^{z_{1}}+e^{z_{2}})=1$.
In this study, specifically, we interpret this probability as the probability
that the first signal is recognized to have a higher frequency by the RNN..
Therefore, ${\rm Softmax}(z_{i})$ represents the probability that the RNN
judge that the $i$-th signal has the higher frequency; If $z_{1}>z_{2}$, the
RNN estimates the first signal to have a higher frequency, and vice versa.

The parameters of this RNN are composed of the weights of the vector, ${\bf
W}^{\rm in}$, and two matrices, ${\bf J}$ and ${\bf W}^{\rm out}$; these
weights are adjusted by training the RNN to solve the task introduced in the
previous subsection. The specific learning procedure is described in the
following subsection. Although, in the reservoir computing schemeMaass  _et
al._ (2002); Sussillo and Abbott (2009) training only changes the output
matrix, ${\bf W}^{\rm out}$, in this study, we adjust the weights of all three
matrices.

###  II.3 Training

To train the RNN for the short-term memory task described above, we adopt a
stochastic gradient descent scheme, which is commonly used in machine-learning
communitiesLeCun  _et al._ (2015). In this scheme, the loss function is first
defined to indicate how far the output of the current RNN is from the correct
answer. This loss function is given as a function of the weights of the
matrices, which provide the parameters of the RNN. The learning process is
carried out by calculating the gradient and optimizing the parameters in the
direction toward which the loss function becomes smaller. Although there are
many optimization algorithms, the basic concept is given by the following
equation:

$W_{t+1}=W_{t}-\eta\nabla_{w}{\mathcal{L}},$ (3)

where $W$ represents the parameters of the RNN (i.e., vector ${\bf W}^{\rm
in}$ and two matrices ${\bf J}$ and ${\bf W}^{\rm out}$), and $\eta$
represents the learning rate. $\nabla_{w}{\mathcal{L}}$ represents the
gradient of the loss function and is calculated via back-propagation through
time (BPTT) Werbos (1990), in which the dynamics of the RNN are first unfolded
in time; then, the derivatives of the loss function are calculated using the
chain-rule. In this study, we adopt the softmax cross-entropy loss function:

${\mathcal{L}}_{\rm CE}=-\sum_{k=1}^{2}\hat{z}_{k}\log z_{k}^{\rm softmax},$
(4)

where $z_{k}^{\rm softmax}\equiv e^{z_{k}}/(e^{z_{1}}+e^{z_{2}})$, and
$\bf\hat{z}$ represents the target label of this task. If $\omega_{1}$ is
larger than $\omega_{2}$, $\hat{\bf z}=(1,0)^{\mathrm{T}}$; otherwise,
$\hat{\bf z}=(0,1)^{\mathrm{T}}$. This target label is defined to satisfy the
demand of this task. As explained, $z_{k}^{\rm softmax}$ gives the probability
that the $k$-th signal is recognized to have higher frequency by the RNN. The
softmax cross-entropy loss, therefore, represents the difference between the
probability estimated by the RNN and the target probability. This loss
function is widely used owing to its computational efficiencyGoodfellow  _et
al._ (2016).

To make the training process stable, L2 norm regularizationBishop (2006) is
applied for the sum of the norm of the RNN weights. Finally, the loss function
is defined as

${\mathcal{L}}={\mathcal{L}}_{\rm CE}+\lambda_{\rm L2}(\sum_{i}W_{{\rm
in},i}^{2}+\sum_{ij}J_{ij}^{2}+\sum_{ij}W_{{\rm out},ij}^{2}),$ (5)

where $\lambda_{\rm L2}$ is set to 0.0001. During the training phase, the loss
function is summed up for 50 input samples; then, the weights of the matrices
are adjusted by the gradient descent. This process is continued for 3,000
iterations. Here we adopt ADAMKingma and Ba (2014) as the specific algorithm
of optimization, which uses $\nabla_{w}{\mathcal{L}}$ calculated by BPTT, as
is mentioned above. The learning rate, $\eta$, is set to 0.001. The training
algorithm is implemented using PyTorchPaszke  _et al._ (2019), which is the
framework for machine learning.

Figure 2:  (a) Fraction judged by the RNN judge to be
$\omega_{2}\geq\omega_{1}$($\circ$). The horizontal axis represents the
difference between $\omega_{1}$ and $\omega_{2}$. Hence, if the fraction is
100% (0%) for $\omega_{2}-\omega_{1}>0(<0)$, respectively, the correct answer
rate ($\times$) is 100%. As shown, if $|\omega_{1}-\omega_{2}|\geq 1$, the
correct answer rate is almost $100\%$, whereas it decreases as
$|\omega_{1}-\omega_{2}|$ decreases. (b) Neural states $\bf x$ at the end of
the second signal ($t=T_{a}$; timing to judge), plotted in the 3-dimensional
principal component (PC) space. The PC axis is computed from the data of the
neural states at $T_{a}$ for various frequency inputs. Colors represent
$\omega_{2}-\omega_{1}$. The plane in the center can separate these neural
states into those satisfying $\omega_{1}>\omega_{2}$.

##  III Results

###  III.1 Short-term memory by transient oscillatory dynamics

After training, 50 pairs of signals having various $\omega_{1}$ and
$\omega_{2}$ were input to the trained RNN. In Fig.2a, the fraction judged by
the RNN to be $\omega_{2}$ was higher than $\omega_{1}$ plotted against
$\omega_{2}-\omega_{1}$. With the condition of $|\omega_{2}-\omega_{1}|>1$,
the accuracy of the choice was greater than 95%. Hence we can see that the RNN
correctly learned to solve the frequency comparison task. To more closely
examine how the neural dynamics compared the frequencies, the neural states
corresponding to various $\omega_{1}$ and $\omega_{2}$ at $t=T_{a}$ (i.e., at
the end of the second signal) were plotted using the three principal
components of $\\{x_{i}\\}$(Fig.2b)Lever  _et al._ (2017). The state changed
continuously in response to the difference in frequencies. Moreover, according
to $\omega_{1}>\omega_{2}$ or $\omega_{1}\leq\omega_{2}$, the states could be
separated by a plane. Hence, the frequency of the first signal was memorized
over the delay period.

Figure 3:  (a) Trajectories of neural activities during a prolonged delay
period ($T_{s}\leq t$, without second signal) were plotted in a 3-dimensional
principal component (PC) space. Trajectories from three different
$\omega_{1}$, giving rise to different states at $t=T_{s}$, are given by
triangle symbol. (b) Neural states $\bf x$ in the PC space at the beginning of
the delay period ($t=T_{s}$) for 50 different $\omega_{1}$’s presented with
different colors. (c) Neural states $\bf x$ in the PC space at the end of the
delay period ($t=T_{f}$) for $\omega_{1}$ corresponding to (b). (d1) Scatter
plot of the norm, $|{\bf x}(t=T_{s})|$, against the frequency of the first
signal, $\omega_{1}$. (d2) Scatter plot of the norm at the end of the delay
period, $|{\bf x}(t=T_{f})|$, against the frequency of the first signal,
$\omega_{1}$. Monotonic dependence is discernible. (e) Neural states $\bf x$
in the PC space at the end of the delay period ($t=T_{f}$) for
$\omega_{1}=1.5,3,4.5$ with various $0\leq\phi\leq\pi$. (f) Directions of
neural trajectories driven by the second signal with $\omega_{2}=3$. The
x-axis shows the amplitude of the trajectory at $t=T_{f}$. Arrows show the
direction and magnitude of the change in the PC3 component of ${\bf
x}(T_{f}+1)-{\bf x}(T_{f})$ for the first input signal with
$\omega_{1}=1.5,2,2.5,3.5,4,4.5$.

Here, the neural activities continued to change over time during the delay
period, as illustrated in 333See Supplemental Material Fig.1 for the
characteristics of the neural activities during the delay period.. It can now
be confirmed that the first signal was memorized by transient dynamics rather
than by the attractor. To further confirm this, the second signal was removed,
and the long-term behavior of neural dynamics after the delay period was
observed. As shown in Fig.3a, when the delay period was prolonged, the neural
states converged to a limit-cycle attractor independent of the signal
frequency, $\omega_{1}$. After the first signal input, neural activities were
attracted to a low-dimensional manifold at which the limit cycle was
located444The fixed-point analysisSussillo and Barak (2013) showed that there
were a few slow points acting as pseudo-saddles through which the neural
activity passed after the first signal. See Supplemental Material Fig.3 for
the slow points plotted in the 3-dimensional PC space.. Then, the neural
activities oscillated with slowly increasing amplitude towards the limit-cycle
attractor. This increase, however, was so slow that the short-term memory was
maintained for a sufficient amount of time. We analyze the duration of memory
in detail below.

Subsequently, we show how memory information (i.e., frequency of signals) was
represented in the neural activity. In Fig.3a, we can see the trend in which
the amplitude of the transient oscillation before the attraction to the limit
cycle has monotonic dependence on $\omega_{1}$. At the end of the delay period
($t=T_{f}$), this trend is remarkable; there is strong negative correlation
between the amplitude, $|{\bf x}(t=T_{f})|$, and the first signal frequency,
$\omega_{1}$(Fig.3d2). Hence, $\omega_{1}$ is encoded by the amplitude of the
transient oscillation(Fig.3c). Notably, this monotonic coding of the input
frequency by the amplitude does not hold immediately after the input of the
first signal. Indeed, at the beginning of the delay period ($t=T_{s}$), there
was no such monotonic dependence on $\omega_{1}$ (see Fig.3b, d1). Therefore
the neural dynamics during the delay period ($T_{s}\leq t\leq T_{f})$ shaped
the manifold from Fig.3b to Fig.3c.

This coding by amplitude is beneficial for discarding information irrelevant
to the task. In the present task, the input signals included phases of the
oscillation apart from the frequency. Hence the RNN should discriminate the
frequency information from the phase information. Here, signals having
different phases were mapped to different points on the same radius trajectory
(Fig.3e). The neural dynamics during the delay period ($T_{s}\leq t\leq
T_{f})$ dampened information on the signal other than its frequency,
$\omega_{1}$.

The neural state corresponding to $\omega_{1}$, encoded by the amplitude of
transient oscillation, was used as the initial state for the response to the
second signal. After the second signal was input, the neural states moved
separately up and down along the PC3-axis according to the sizes of
$\omega_{2}$ and $\omega_{1}$. If $\omega_{2}>\omega_{1}$, it moved in the
positive direction along the PC3-axis, and vice versa. In Fig.3f, changes in
the PC3 component of ${\bf x}(t)$ after the second signal with $\omega_{2}=3$
are plotted against the amplitude of ${\bf x}(t)$ oscillation at $t=T_{f}$. As
shown, depending on the amplitude of the oscillation at $t=T_{f}$, PC3 moved
upwards when $|{\bf x}(t=T_{f})|$ was small (corresponding to
$\omega_{2}<\omega_{1}$), and it moved downwards when $|{\bf x}(t=T_{f})|$ was
large (corresponding to $\omega_{2}>\omega_{1}$). This property separated the
neural states in the direction of PC3, as shown in Fig.2b, and it enabled the
RNN to correctly solve the frequency comparison task by using the short-term
memory for $\omega_{1}$, as coded in the transient amplitude.

To verify the generality of short-term memory encoded by the amplitude of
transient oscillation, we examined two other settings. First, we changed the
length of the delay period as follows. During the training phase, $T_{d}$ was
homogeneously distributed as $T_{d}\in[75,105]$, and during the test phase,
$T_{d}$ was fixed at 90. With this longer delay period setting, we confirmed
that the present mechanism also worked 555See Supplemental Material Fig.2 for
the neural activity of trained RNN with longer delay period.. Second, we
trained the RNN for two other comparison tasks, which requested we compare the
velocity and noise variance. In the former task, the first and second signals
were given by $u_{1,2}(t)=a_{1,2}t+b+\eta_{1,2}(t)$, and the RNN was required
to determine which of the velocities, $a_{1,2}$, was larger. In the latter
task, the first and second signals were given by
$u_{1,2}=\sin(t+\phi)+\eta_{1,2}(t)$, where $\eta_{1,2}(t)$ was a random
Gaussian variable having an average of zero and a standard deviation of
$\sigma_{\rm ex}^{1,2}$. The RNN was then required to determine which of the
variances, $\sigma_{\rm ex}^{1,2}$, was larger. For both the tasks, a
monotonic dependence between the L2 norm of the neural states (the amplitude
of transient oscillation) and the parameter to be compared was revealed at the
end of the delay period ($t=T_{f}$)666See Supplemental Material Fig.4 and 5
for the results of the different tasks.. The information of the velocity,
$a_{1}$, (for the former) and the noise variance, $\sigma_{\rm ex}^{1}$, (for
the latter) was memorized as the amplitude of transient oscillation.

###  III.2 Convergence to the limit cycle

As described previously, transient oscillatory dynamics eventually converged
to the limit cycle. When this trajectory converged, the memory was forgotten,
and the information of the input was lost. Hence, for this memory to work and
to accomplish the task, the time to converge to the limit cycle must be
sufficiently long. As a measure of convergence of ${\bf x}(t)$ to the limit
cycle, we first introduced $\ell(t)=|{\bf x}(t)-{\bf x}(t+T_{L})|$, where
$T_{L}$ is the period of the limit cycle. Subsequently, we estimated the
convergence time to the limit cycle as the time, $t$, when $\ell(t)$ was
sufficiently small (i.e., the time $t$ at which $\ell(t)$ first satisfies
$\ell(t)\leq 0.05$ for the first time). We computed the average convergence
time over 100 trajectories, ${\bf x}(t)$, for different $\omega_{1}$ values
and noted that the convergence time was a dozen times longer than $T_{f}$. We
also found that during training, the convergence time increased alongside an
increase in the accuracy of the task 777See Supplemental Material Fig.6 for
the convergence time to the limit cycle.. These results suggest that the
short-term memory was maintained over a sufficiently long period, prolonged
after the delay period, and the length of the convergence time was related to
the performance of informational processing.

###  III.3 Robustness to noise

As described in the introduction, the mechanism of the robustness to noise in
the short-term memory by transient oscillation remains elusive. Hence, the
robustness of the memory was examined by adding the noise term in Eq.(1) using
the Langevin equation:

$\dot{x}_{i}=-x_{i}+\sum_{j=1}^{N}J_{ij}\tanh(x_{j})+W_{i}^{\rm in}u+\xi_{i},$
(6)

where $\xi_{i}$ is a random Gaussian variable with an average of zero and a
standard deviation of $\sigma_{\rm neu}$. As with the no-noise setting, we
discretized the dynamics and obtained the following equation:
$x_{i}(t+1)=(1-\alpha)x_{i}(t)+\alpha(\sum_{j=1}^{N}J_{ij}\tanh(x_{j}(t))+W_{i}^{in}u(t))+\sqrt{\alpha}\xi_{i}(t)$.
The RNN was trained under a noise level, $\sigma_{\rm neu}^{\rm train}$, with
regularization for the average squared norm of neural activity,
$\frac{1}{T_{a}}\sum_{t}^{T_{a}}\sum_{i}^{N}x_{i}(t)^{2}$. Regularization was
applied to the loss function in the form of

$\mathcal{L}_{\rm noise}={\mathcal{L}}+\lambda_{\rm
act}\frac{1}{T_{a}}\sum_{t}^{T_{a}}\sum_{i}^{N}x_{i}(t)^{2}.$ (7)

The reason that we introduced this regularization term is as follows. Notably,
if the average squared norm,
$\frac{1}{T_{a}}\sum_{t}^{T_{a}}\sum_{i}^{N}x_{i}(t)^{2}$, increases,
robustness will increase because the derivative of $\tanh(x)$ will generally
decrease. By using this regularization term to suppress the norm of neural
activity, we can avoid such trivial robustness. We experimentally determined
$\lambda_{\rm act}$ so that the norm of the internal dynamics of RNNs trained
with noise ($\sigma_{\rm neu}^{\rm train}=0.04$) would be comparable to those
trained without noise ($\sigma_{\rm neu}^{\rm train}=0$), and we adopted
$\lambda_{\rm act}=30$. Subsequently, the trained network was tested for
solving the task under the noise level, $\sigma_{\rm neu}^{\rm test}$.

Figure 4:  (a) Accuracy of the frequency comparison task plotted against the
noise level, $\sigma_{\rm neu}^{\rm test}$. Different colors show the RNN
trained at different levels of noise ($\sigma_{\rm neu}^{\rm train}=0,0.04)$.
Bold lines show the mean score, and shaded areas show the standard deviations.
(b) The cumulative contribution ratio of the perturbed trajectories up to the
10th principal component (i.e. the cumulative percentage of the eigenvalues
corresponding to 1st-to-10th eigenvectors
$(\sum_{i=1}^{10}\lambda_{i})/(\sum_{i=1}^{256}\lambda_{i})$. These are
plotted against the time following the application of perturbation for
$\sigma_{\rm neu}^{\rm train}=0$(blue), and $\sigma_{\rm neu}^{\rm
train}=0.04$(red). If it is larger, more points are restricted within the
lower-dimensional manifold. Because the perturbation is completely random, the
perturbation dimension is high immediately after it is applied. However, over
time, it is attracted to a lower-dimensional manifold, especially for
$\sigma_{\rm neu}^{\rm train}=0.04$.

As depicted in Fig.4a, the accuracy of the short-term memory task decreased
with the applied noise, $\sigma_{\rm neu}^{\rm test}$, for the network trained
without noise. In contrast, for the model trained with a sufficient noise
level, $\sigma_{\rm neu}^{\rm train}$, the drop in the accuracy by the
increase in noise level was suppressed, even up to a noise level higher than
$\sigma_{\rm neu}^{\rm train}$. The system gained a higher robustness to
noise, which was beyond the level added during learning. Notably, the RNN
trained with neural noise realized the same mechanism of the short-term memory
with transient oscillation.

To examine the achievement of the stability of short-term memory, the
following perturbation analysis was performed. We set $\sigma_{\rm neu}^{\rm
test}=0$, and subsequently, we introduced an instantaneous perturbation during
the delay period (i.e., $T_{s}\leq t_{\rm per}\leq T_{f}$, as ${\bf x}_{\rm
per}(t)={\bf x}(t)+\bf\delta x\delta_{t,t_{\rm per}}$), where ${\bf\delta x}$
is a random variable that follows a Gaussian distribution of mean zero and of
variance $\sigma_{\rm per}^{2}$, and $\delta_{t,t_{\rm per}}$ represents a
Kronecker delta.

Intuitively, one might expect that the distance between the perturbed
trajectory and the original trajectory would be decreased for the networks
trained under noise. Here, we first computed the L2 norm distance between the
perturbed and original trajectories: $L_{per}=|{\bf x}_{per}(t)-{\bf x}(t)|$.
In contrast to the naive expectation, however, the distance does not depend
substantially on the trained noise level, $\sigma_{\rm neu}^{\rm train}$.
Hence, the distance between the trajectories cannot explain the memory
robustness for a system trained with noise $\sigma_{\rm neu}^{\rm train}$. To
understand the robustness, the direction of the perturbation of the
trajectories must be considered. Indeed, the shift of the trajectory along the
original trajectory was not harmful to the memory discussed, because the
amplitude of the transient oscillation was not affected by the shift. Hence,
the dimension of the manifold spanned by the perturbed trajectories is more
important to the robustness of memory than the distances between the perturbed
trajectories and the original trajectories.

Figure 5:  Histogram of the frequency distribution of the real part of the
eigenvalues of $\exp[\frac{1}{30}\sum_{t=15}^{45}{\bf G}(t)]$. We calculated
${\bf G}(t)$ for 50 samples of the input signals. Notably, there were 256
eigenvalues for one $\exp[\frac{1}{30}\sum_{t=15}^{45}{\bf G}(t)]$. Blue is
the result for the RNN trained without noise, and red is the result for the
RNN trained with noise($\sigma_{\rm neu}^{\rm train}=0.04$). In the latter,
the eigenvalues were biased to the negative side overall.

Accordingly, we estimated the dimension of the manifold spanned by a large
number of perturbed trajectories. We adopted a principal component analysis
(PCA) for an ensemble of these perturbed trajectories. If perturbation causes
a shift mainly along the trajectory direction, perturbed trajectories should
be restricted within a low-dimensional manifold along the original trajectory.
Hence, the dimension of perturbed trajectories estimated by PCA will work as a
measure of robustness in the directions of perturbed trajectories. To estimate
the dimension, we computed the percentage of the variance corresponding the
1st-to-10th eigenvectors for each time after the perturbation as
$(\sum_{i=1}^{10}\lambda_{i})/(\sum_{i=1}^{256}\lambda_{i})$, where
$\lambda_{i}$ represents the eigenvalue of covariance matrix with an ordering
to satisfy $\lambda_{1}\geq\lambda_{2}\geq...\geq\lambda_{256}$.

This percentage was larger when the perturbed trajectory was restricted to a
lower-dimensional manifold. As shown in Fig. 4, for the RNN trained with noise
$\sigma_{neu}$, the trajectory falls into a lower-dimensional manifold with
over time following the perturbation. The percentage for such robust RNN is
larger than the RNN trained without neural noise. This result implies that the
dynamics after the perturbation were more restricted to a lower-dimensional
manifold for the model trained with noise. This low-dimensional compression is
key to the robustness of the neural noise. Notably, another method of
estimating the attractor dimension was proposed in Sirovich (1989); Ciliberto
and Nicolaenko (1991), where the number of eigenvectors of the principal
components required to achieve 90% of the total variance was calculated. We
have confirmed that similar results have been obtained by applying this method
888See Supplemental Material Fig.7 for the estimated dimension of perturbed
trajectories by the method in Sirovich (1989); Ciliberto and Nicolaenko
(1991).

To characterize the mechanism of this convergence of the transient trajectory
to low-dimensional space, we focused on the Jacobi matrix, ${\bf G}(t)$, of
this system. The Jacobi matrix is defined as the following equation:

$G_{ij}(t)=J_{ij}\tanh^{\prime}(x_{j}(t)).$ (8)

The perturbed dynamics convergence can be estimated by the eigenvalues of
matrix $e^{\frac{1}{T}\int dt{\bf G}(t)}$ along the trajectory, as was adopted
in the calculation of finite-time Lyapunov exponentsCrisanti  _et al._ (1988).
Therefore, we investigated the difference of eigenvalues between the RNNs
trained with $\sigma_{\rm neu}^{\rm train}=0$ and those trained with
$\sigma_{\rm neu}^{\rm train}=0.04$. Because the value of the Jacobi matrix
depends on the specific trajectory, ${\bf x}(t)$, we calculated the histograms
of the eigenvalues for 50 samples of input signals (Fig.5). We also adopted
the discretization to calculate $e^{\frac{1}{T}\int dt{\bf G}(t)}$ as
$\exp[\frac{1}{30}\sum_{t=15}^{45}{\bf G}(t)]$.

For the latter RNN, there were more eigenvalues with negative real parts. As
shown in Fig.5, the number of eigenvalues with positive real parts was
reduced. This result suggests that the restriction to a low-dimensional
manifold in the robust RNN was caused by the compression of many modes of
perturbed dynamics. This robustness caused by the contraction of transient
trajectories was clearly distinguishable from that caused by the stability of
fixed-point attractors.

###  III.4 Different types of neural networks

To investigate the range of validity of transient oscillatory short-term
memory, we considered two different types of neural networks. First, we
considered neural networks in which synapses obey Dale’s law (i.e., the neural
network consists of excitatory and inhibitory neurons, and synapses extended
from a given neuron are either all excitatory , or inhibitory)Dale (1935);
Song  _et al._ (2016). Second, we considered neural networks in which synapses
were sparse. We call these networks excitatory-inhibitory and Sparse networks,
respectively.

We trained the RNN to solve the frequency comparison task while satisfying
each condition. Under the excitatory-inhibitory networks condition, as in
Orhan and Ma (2017), hidden neurons were divided into excitatory and
inhibitory neurons in a 4-to-1 ratio. $\bf J$ was constrained so it can be
decomposed as ${\bf J}={\bf J}^{+}{\bf D}$, where ${\bf J}^{+}$ is a matrix
with all components greater than or equal to zero, and $\bf D$ is a diagonal
matrix that satisfies $D_{ii}=1$ if $i\leq 180$ and $D_{ii}=-1$ if $i>180$.
With this constraint, it can be said that $\bf J$ obeys Dale’s law. During the
training phase, ${\bf D}$ was fixed, and ${\bf W}^{\rm in}$, ${\bf J}^{+}$ and
${\bf W}^{\rm out}$ were adjusted. In Sparse networks, we set the percentage
of synapses with nonzero weight to approximately 20%. We adopted the deep-
rewiring techniqueBellec  _et al._ (2017) to train under constraints where
only 20% of synapses have nonzero values, and all others are zero. As in
Supplementary Material Figs. 8 and 9, after learning, the RNNs solved the
frequency comparison task by maintaining short-term memory by transient
oscillatory dynamics 999See Supplemental Material Fig.8 and 9 for the results
under the excitatory-inhibitory networks condition and the sparse networks
condition.

##  IV Discussion

In this study, we uncovered a generic scheme for short-term memory sustained
by transient oscillatory dynamics by training RNNs to achieve a task requiring
short-term memory to compare two interspaced input signalsBarak (2017);
Richards  _et al._ (2019). We demonstrated that short-term memory was encoded
in the amplitude of transient oscillations of neural activities. The neural
state given by high-dimensional dynamical systems fell into a low-dimensional
manifoldSussillo and Barak (2013) and exhibited transient oscillation, which
slowly approached a limit-cycle attractor. With the passage of time,
continuous information from the first signal input, such as the frequency and
velocity of the input signal, was encoded in the amplitude of the transient
oscillation and maintained as short-term memory. Other irrelevant information
in the inputs (e.g., phase and noise) was discarded during the transient
dynamics. Hence, short-term memory is encoded and maintained by transient
dynamics and is robust to external noise.

The proposed mechanism, wherein the memory is encoded in transient
oscillation, contrasts with the view of memory as fixed-point attractors
(i.e., memories as attractors). When the RNN is trained to store discrete
information (e.g., the possible input signal frequency candidates are limited
to 1, 2, 3, 4, or 5 Hz), the memories are encoded by multiple fixed-point
attractors (i.e., persistent activities)Christian  _et al._ (2005). In the
task adopted in this study, by contrast, storing continuous information was
needed. Hence, encoding into the amplitude of transient oscillation occurred.
Notably, as an alternative mechanism to encoding continuous information, a
line attractor was proposedSeung (1996); Mante and et al. (2013). Here,
attractors coexisted continuously on a line in the state space of neural
activities, along which the fixed-point was marginally stable (i.e., one of
the eigenvalues of the Jacobi matrix (with the eigenvector along the line)
should be zero). However, in autonomous dynamical systems, such marginally
stable attractors are not generic, and for their existence, special
constraints are required. Indeed, in this study, such line attractors were not
shaped by learning101010For the possible applicability of chaotic attractors
to short-term memory, see J. S. Nicolis and I. Tsuda, Chaotic dynamics of
information processing: The “magic number seven plus-minus two” revisited,
Bulletin of Mathematical Biology47, 343(1985)..

As another mechanism of short-term memory, synaptic plasticity has been
proposedMongillo  _et al._ (2008); Mi  _et al._ (2017); Taher  _et al._
(2020). In this case, the synaptic connections are modeled as time-varying not
only during training, but also during inferencing. The memory information is
encoded in the time-varying strength of the synapses. In this study, because
we focused on neural-based short-term memory, synaptic connections (i.e., the
vector ${\bf W}^{\rm in}$ and two matrices ${\bf J}$ and ${\bf W}^{\rm out}$)
change only during training and were fixed thereafter. In future research, we
plan to investigate how short-term memory based on transient oscillation can
be realized when considering synaptic plasticity after training.

It remains to be seen if the present scheme for short-term memory can be
adopted biologically. The results suggest that the short-term memory encoded
by the amplitude of transient oscillation works over a wide range of systems.
This may be plausible, because the present mechanism is generally represented
in terms of dynamical systems, and it is not difficult to generate oscillatory
dynamics from an ensemble of neurons. Here, we offer two remarks. We expect
that the present scheme is also valid therein, as it is robust to noise in
neural dynamics. Next, oscillatory dynamics from certain modes of neural
activities evoked by inputs were observed from neural data referencesMuller
_et al._ (2018); H.  _et al._ (2018). Notably, such dynamics are not
necessarily observed just by computing the average neural activity, which is
often a rather stationary reference independent of inputs. Indeed, this is
true in the proposed model. The prominent oscillation depending on the input,
as shown in Fig.2, is observed only by taking appropriate principal
components. If we compute just the average activity over all neurons, it is
almost stationary, and the oscillation is difficult to discern, which is
consistent with experimental observations.

In this study, attraction to the low-dimensional manifold of perturbation
allowed for the transient trajectory to be robust to noise. Recently, dynamic
robustness has been regarded as an important property of information
processing in the brainWarasinee  _et al._ (2017). The present results provide
general insight into the mechanism of this type of robustness. Furthermore, in
cell and developmental systems, such dynamic robustness has also been
discussed as homeorhesisWaddington (1957); Matsushita and Kaneko (2020); Young
_et al._ (2017), and the results may be applicable to such phenomena.

###### Acknowledgements.

The authors would like to express their gratitude to Tomoki Kurikawa, Ichiro
Tsuda, Hiromichi Suetani, and Tetsuhiro S. Hatakeyama for their comments and
discussions. This research was partially supported by a Grant-in-Aid for
Scientific Research (A) (20H00123) from the Japanese Society for the Promotion
of Science (JSPS).

## References

  * Jonides  _et al._ (2008) J. Jonides, R. L. Lewis, and et al., Annual Review of Psychology 59, 193 (2008).
  * Colby  _et al._ (1996) C. L. Colby, J. R. Duhamel, and M. E. Goldberg, Journal of Neurophysiology 76, 2841 (1996), pMID: 8930237, https://doi.org/10.1152/jn.1996.76.5.2841 .
  * Funahashi  _et al._ (1989) S. Funahashi, C. Bruce, and P. Goldman-Rakic, Journal of neurophysiology 61, 331 (1989).
  * Fuster and Alexander (1971) J. Fuster and G. Alexander, Science 173, 652 (1971).
  * Wang (2001) X. J. Wang, Trends in Neurosciences 24, 455 (2001).
  * Haider and McCormick (2009) B. Haider and D. A. McCormick, Neuron 62, 171 (2009).
  * Baeg  _et al._ (2003) E. H. Baeg, Y. B. Kim, K. Huh, I. Mook-Jung, H. T. Kim, and M. W. Jung,  _Neuron_ , Neuron 40, 177 (2003).
  * Bondanelli and Ostojic (2020) G. Bondanelli and S. Ostojic, PLoS Comput Biol 16 (2020).
  * Druckmann and Chklovskii (2012) S. Druckmann and D. Chklovskii, Current biology : CB 22 (2012), 10.1016/j.cub.2012.08.058.
  * Goldman (2009) M. S. Goldman, Neuron 61, 621 (2009).
  * Orhan and Pitkow (2020) E. Orhan and X. Pitkow, in  _International Conference on Learning Representations_ (2020).
  * Orhan and Ma (2019) E. Orhan and W. J. Ma, Nature Neuroscience 22, 275 (2019).
  * Rabinovich  _et al._ (2008) M. Rabinovich, R. Huerta, and G. Laurent, Science 321, 48 (2008).
  * Sussillo and Barak (2013) D. Sussillo and O. Barak, Naural Comput. 25(3), 626 (2013).
  * Maheswaranathan  _et al._ (2019) N. Maheswaranathan, A. Williams, and et al, in  _Advances in Neural Information Processing Systems 32_ (Neural Information Processing Systems Foundation, 2019) pp. 15629–15641.
  * Romo  _et al._ (1999) R. Romo, C. Brody, A. Hernández, and L. L., Nature 399, 470 (1999).
  * Rumelhart  _et al._ (1986) D. E. Rumelhart, G. E. Hinton, and R. J. Williams, Nature 323, 533 (1986).
  * Werbos (1990) P. J. Werbos, Proceedings of the IEEE 78, 1550 (1990).
  * Rivkind and Barak (2017) A. Rivkind and O. Barak, Phys. Rev. Lett. 118, 258101 (2017).
  * Note (1) Because $x_{i}$ takes both positive and negative values, it does not represent the neural activity per se. Instead, in the context of neuroscience, $(1+\mathop{tanh}\nolimits(x_{i}))/2$ is often regarded as the firing rate of the neuron.
  * Bishop (2006) C. M. Bishop,  _Pattern Recognition and Machine Learning_ (Springer, 2006).
  * Note (2) The relationship between logit $z_{1}-z_{2}$ and the probability calculated by the softmax function is as follows. By definition of the softmax function, we have $p_{1}=e^{z_{1}}/(e^{z_{1}}+e^{z_{2}})=1/(1+e^{-(z_{1}-z_{2})})$. Then, by the definition of logit function, we get the logit of $p_{1}$ as $\mathop{log}\nolimits(p_{1}/(1-p_{1}))=\mathop{log}\nolimits(e^{z_{1}-z_{2}})=z_{1}-z_{2}$. From the above explanation, the output of the softmax function can be treated as a probability because $p_{1}$ and $p_{2}$ satisfy $0\leq p_{1},p_{2}\leq 1$, and $p_{1}+p_{2}=e^{z_{1}}/(e^{z_{1}}+e^{z_{2}})+e^{z_{2}}/(e^{z_{1}}+e^{z_{2}})=1$. In this study, specifically, we interpret this probability as the probability that the first signal is recognized to have a higher frequency by the RNN.
  * Maass  _et al._ (2002) W. Maass, T. Natschläger, and H. Markram, Neural Comput. 14(11) (2002).
  * Sussillo and Abbott (2009) D. Sussillo and L. Abbott, Neuron 63, 544 (2009).
  * LeCun  _et al._ (2015) Y. LeCun, Y. Bengio, and G. Hinton, Nature 521, 436 (2015).
  * Goodfellow  _et al._ (2016) I. Goodfellow, Y. Bengio, and A. Courville,  _Deep Learning_ (MIT Press, 2016) http://www.deeplearningbook.org.
  * Kingma and Ba (2014) D. Kingma and J. Ba, International Conference on Learning Representations (2014).
  * Paszke  _et al._ (2019) A. Paszke, S. Gross, and et al, in  _Advances in Neural Information Processing Systems 32_ (Neural Information Processing Systems Foundation, 2019) pp. 8024–8035.
  * Lever  _et al._ (2017) J. Lever, M. Krzywinski, and N. Altman, Nat Methods 14, 641 (2017).
  * Note (3) See Supplemental Material Fig.1 for the characteristics of the neural activities during the delay period.
  * Note (4) The fixed-point analysisSussillo and Barak (2013) showed that there were a few slow points acting as pseudo-saddles through which the neural activity passed after the first signal. See Supplemental Material Fig.3 for the slow points plotted in the 3-dimensional PC space.
  * Note (5) See Supplemental Material Fig.2 for the neural activity of trained RNN with longer delay period.
  * Note (6) See Supplemental Material Fig.4 and 5 for the results of the different tasks.
  * Note (7) See Supplemental Material Fig.6 for the convergence time to the limit cycle.
  * Sirovich (1989) L. Sirovich, Physica D: Nonlinear Phenomena 37, 126 (1989).
  * Ciliberto and Nicolaenko (1991) S. Ciliberto and B. Nicolaenko, Europhysics Letters (EPL) 14, 303 (1991).
  * Note (8) See Supplemental Material Fig.7 for the estimated dimension of perturbed trajectories by the method in Sirovich (1989); Ciliberto and Nicolaenko (1991).
  * Crisanti  _et al._ (1988) A. Crisanti, G. Paladin, and A. Vulpiani, J Stat Phys 53, 583 (1988).
  * Dale (1935) H. Dale, Proc R Soc Med. 28(3), 319 (1935).
  * Song  _et al._ (2016) H. F. Song, G. R. Yang, and X.-J. Wang, PLOS Computational Biology 12, 1 (2016).
  * Orhan and Ma (2017) A. Orhan and W. Ma, Nat Commun 8 (2017).
  * Bellec  _et al._ (2017) G. Bellec, D. Kappel, W. Maass, and R. A. Legenstein, CoRR abs/1711.05136 (2017).
  * Note (9) See Supplemental Material Fig.8 and 9 for the results under the excitatory-inhibitory networks condition and the sparse networks condition.
  * Barak (2017) O. Barak, Curr Opin Neurobiol. 46, 1 (2017).
  * Richards  _et al._ (2019) B. Richards, T. Lillicrap, P. Beaudoin, and et al., Nat Neurosci 22, 1761 (2019).
  * Christian  _et al._ (2005) K. M. Christian, R. Ranulfo, and D. B. Carlos, science 307, 1121 (2005).
  * Seung (1996) H. S. Seung, Proceedings of the National Academy of Sciences 93, 13339 (1996).
  * Mante and et al. (2013) V. Mante and et al., Nature 503, 78 (2013).
  * Note (10) For the possible applicability of chaotic attractors to short-term memory, see J. S. Nicolis and I. Tsuda, Chaotic dynamics of information processing: The “magic number seven plus-minus two” revisited, Bulletin of Mathematical Biology47, 343(1985).
  * Mongillo  _et al._ (2008) G. Mongillo, O. Barak, and M. Tsodyks, Science 319, 1543 (2008).
  * Mi  _et al._ (2017) Y. Mi, M. Katkov, and M. Tsodyks, Neuron 93, 323—330 (2017).
  * Taher  _et al._ (2020) H. Taher, A. Torcini, and S. Olmi, PLOS Computational Biology 16, 1 (2020).
  * Muller  _et al._ (2018) L. Muller, F. Chavane, J. Reynolds, and et al., Nat Rev Neurosci 19, 255 (2018).
  * H.  _et al._ (2018) Z. H., W. AJ., and et al., Neuron 27, 1269 (2018).
  * Warasinee  _et al._ (2017) C. Warasinee, W. Xiao-Jing, and et al., Neuron 93, 1504 (2017).
  * Waddington (1957) C. H. Waddington,  _The strategy of the genes: A discussion of some aspects of theoretical biology_ (Allen & Unwin, London, 1957).
  * Matsushita and Kaneko (2020) Y. Matsushita and K. Kaneko, Phys. Rev. Research 2, 023083 (2020).
  * Young  _et al._ (2017) J. T. Young, T. S. Hatakeyama, and K. Kaneko, PLOS Computational Biology 13, 1 (2017).

Supplemental Material

Figure S-1: (a)The neural activity of one neuron $x_{100}(t)$ during the delay
period. Different colors represent the activity for different $\omega_{1}$.
(b) The firing rate of neurons during the delay period(Raster plot). Neurons
are sorted by the time of their peak firing. One can see the propagation of
neural activities, whose ’speed’ slows down during the delay period. As shown
in the figure, the neural activity changes dynamically during the delay
period. This means that stable short-term memories are encoded by transient
neural dynamics. The dynamics slow down during the delay period so that the
short-term memory is maintained for a sufficient time. Note that although
neural activities slow down, they converge not to a fixed point, but to the
limit cycle.  Figure S-2: The neural activity of trained RNN with the
condition that in the training phase, $T_{d}$ is homogeneously distributed as
$T_{d}\in[75,105]$ and in the test phase, $T_{d}$ is fixed at 90.(c.f. in main
paper, in the training phase, $T_{d}$ is homogeneously distributed as
$T_{d}\in[25,35]$ and in the test phase, $T_{d}$ is fixed at 30).
(a)Trajectories of neural activities during prolonged delay period($T_{s}\leq
t$, without second signal), were plotted in 3-dimensional PC space.
Trajectories from three different $\omega_{1}$, giving rise to different
states at $t=T_{s}$, were given by triangle symbol.(b) Neural states $\bf x$
in the PC space at the beginning of the delay period($t=T_{s}$) for 50
different $\omega_{1}$’s presented in different colors. (c)Neural states $\bf
x$ in the PC space at the end of the delay period($t=T_{f}$) for $\omega_{1}$
corresponding to (b). (d1) Scatter plot of the norm $|{\bf x}(t=T_{s})|$
against the frequency of the first signal $\omega_{1}$.(d2) Scatter plot of
the norm at the end of the delay period $|{\bf x}(t=T_{f})|$ against the
frequency of the first signal $\omega_{1}$. Monotonic dependence is
discernible. (e)Neural states $\bf x$ in the PC space at the end of the delay
period($t=T_{f}$) for $\omega_{1}=1.5,3,4.5$. with various $0\leq\phi\leq\pi$.
As in the case of $T_{d}=30$(main paper, Fig.3), a transient oscillation is
formed during the delay period, and at $t=T_{f}$, $\omega_{1}$ is encoded as
amplitude of transient oscillation ($|{\bf x}(t=T_{f})|$ has a monotonic
dependence on $\omega_{1}$).  Figure S-3: ’Slow points’ plotted in the
3-dimensional PC space with eigenvectors corresponding to unstable modes. The
speed at a point on the trajectory $|{\bf\dot{x}}|_{2}=|-{\bf x}+{\bf
J}\tanh({\bf x})|_{2}$ was computed, and its local minimum points were
obtained using the gradient method. The fixed points where the speed is
completely zero are only the central points of the limit cycle, but there are
some ’slow points’ where the speed is very less, as depicted in the figure.
These ’slow points’ are also useful to understand the nature of the dynamics
because they often work as pseudo stable plane or pseudo saddles; the
eigenvalues of the Jacobian around these ’slow points’ were calculated. It was
found that the trajectory is attracted to these ’slow points’, because most of
the modes are negative(only one mode is positive, so they work as pseudo
saddles). Indeed, the timing of neural activity slows down(Fig.S-7b) coincides
with the timing of passage through these ’slow points’, indicating that they
contribute to the formation of the slow manifold.  Figure S-4:  RNN was
trained for velocity comparison task in which the first and second signals are
given by $u_{1,2}(t)=a_{1,2}t+b+\eta$. Then, it was asked to determine which
of the velocities $a_{1,2}$ is larger. $\eta$ is the Gaussian white noise.
(a)Sample signals. Different colors show different values of $a_{1}$.
(b)Neural activity of trained RNN plotted in 3-dimensional PC space. Color
points show the neural state at the end of the delay period ($t=T_{f}$).
Oscillatory neural activity, which converged to a limit cycle after a long
period of time, was observed during the delay period, (c) Scatter plot of the
norm at the end of the delay period relative to the base state $|{\bf x}-{\bf
x}_{0}|$ plotted against the frequency of the first signal $a_{1}$. The base
state corresponds to the neural activity for the input signal’s velocity which
is $a_{1}=-0.15$. The input signal information is embedded in the amplitude of
the transient oscillation.  Figure S-5: RNN was trained for noise variance
comparison task. Here, the first and second signals are given by
$u_{1,2}=\sin(t+\phi)+\eta$ in which $\eta$ as random Gaussian variable
following $\eta\sim\mathcal{N}(0,\ \sigma_{1,2}^{in})$.Then, it is asked which
of the variance $\sigma_{1,2}$ is larger.
(a)Sample signals. Different colors show different $\sigma_{1}^{in}$.
(b)Neural activity of trained RNN plotted in 3-dimensional PC space. Color
points show the neural state at the end of the delay period ($t=T_{f}$).
Oscillatory neural activity was observed during the delay period, which
converged to a limit cycle after a long period of time. (c) Scatter plot of
the norm at the end of the delay period relative to the base state $|{\bf
x}-{\bf x}_{0}|$ plotted against the frequency of the first signal $a_{1}$.
The base state corresponds to the neural activity for the input signal
velocity with $\sigma_{in}=0$. The input signal information is embedded in the
amplitude of the transient oscillation. Figure S-6: We estimated the
convergence time to the limit cycle as the time $t$ when $|{\bf x}(t)-{\bf
x}(t+T_{L})|$ was sufficiently small.
(a)$|{\bf x}(t)-{\bf x}(t+T_{L})|$ with the time from $t=T_{s}$. Each of the
three lines corresponds to the results of the RNN model under training
(Epoch=1000, 2000, 3000). Since the delay period is 30, the convergence time
is sufficiently larger than the delay period in all cases. In addition, the
convergence time increases as training progresses. (b) Time course of
learning. It reaches sufficient accuracy before Epoch 1000, and from there on,
the accuracy slowly improves as learning progresses.  Figure S-7: We estimated
dimension of perturbed trajectories by the method in [35, 36], i.e. the number
of eigenvectors of principal components need to give 90% of the total
variance. This value is smaller when the perturbed trajectory is restricted to
a lower-dimensional manifold. The number of eigenvectors is plotted against
the time after the application of perturbation, for $\sigma_{\rm neu}^{\rm
train}=0$(blue), and $\sigma_{\rm neu}^{\rm train}=0.04$(red). As in the
result of main paper Fig.4., the dynamics after the perturbation are more
restricted to a lower-dimensional manifold for the model trained with noise.
Figure S-8: We trained RNN to solve the frequency comparison task while
satisfying E-I networks condition. In this condition, hidden neurons were
divided into excitatory-inhibitory neurons in a 4-to-1 ratio. Excitatory
neurons only have ’excitatory synapes’ which has positive value of synaptic
weight, and Inhibitory neurons only have ’inhibitory synapes’ which has
negative value of synaptic weight. As in the result of non constrained
RNN(main paper), there are transient oscillation and amplitude coding of input
signals.
(a) Scatter plot of the norm at the end of the delay period $|{\bf
x}(t=T_{f})|$ against the frequency of the first signal $\omega_{1}$.
Monotonic dependence is discernible. (b) Neural states $\bf x$ in the PC space
at the end of the delay period($t=T_{f}$) for different $\omega_{1}$.  Figure
S-9: We trained RNN to solve the frequency comparison task while satisfying
Sparse networks condition. In this condition, the percentage of synapses with
non-zero weight is set to 20%. As in the result of non constrained RNN(main
paper), there are transient oscillation and amplitude coding of input signals.
(a) Scatter plot of the norm at the end of the delay period $|{\bf
x}(t=T_{f})|$ against the frequency of the first signal $\omega_{1}$.
Monotonic dependence is discernible. (b) Neural states $\bf x$ in the PC space
at the end of the delay period($t=T_{f}$) for different $\omega_{1}$.