# Multiphysics Numerical Method for Modeling Josephson Traveling-Wave
Parametric Amplifiers

Samuel T. Elkin, Michael Haider, and Thomas E. Roth Manuscript received XX XX,
2024; revised XX XX, 2024.Samuel T. Elkin and Thomas E. Roth are with the
Elmore Family School of Electrical and Computer Engineering, Purdue
University, West Lafayette, IN 47907 USA and the Purdue Quantum Science and
Engineering Institute, West Lafayette, IN 47907 USA (e-mail:
rothte@purdue.edu)Michael Haider is with the TUM School of Computation,
Information and Technology, Technical University of Munich, 85748 Garching,
Germany

###### Abstract

Josephson traveling-wave parametric amplifiers (JTWPAs) are wideband,
ultralow-noise amplifiers used to enable the readout of superconducting
qubits. While individual JTWPAs have achieved high performance, behavior
between devices is inconsistent due to wide manufacturing tolerances.
Amplifier designs could be modified to improve resilience towards variations
in amplifier components; however, existing device models often rely on
analytical techniques that typically fail to incorporate component variations.
To begin addressing this issue, a 1D numerical method for modeling JTWPAs is
introduced in this work. The method treats the Josephson junctions and
transmission lines in an amplifier as coupled subsystems and can easily
incorporate arbitrary parameter variations. We discretize the transmission
line subsystem with a finite element time domain method and the Josephson
junction subsystem with a finite difference method, with leapfrog time
marching used to evolve the system in time. We validate our method by
comparing the computed gain to an analytical model for a traditional JTWPA
architecture and one with resonant phase matching. We then use our method to
demonstrate the impact of variations in Josephson junctions and phase-matching
resonators on amplification. In future work, the method will be adjusted to
incorporate additional amplifier architectures and extended to a 3D full-wave
approach.

###### Index Terms:

Multiphysics modeling, computational electromagnetics, Josephson junction,
Josephson traveling-wave parametric amplifier, resonant phase matching.

##  I Introduction

Superconducting qubits are a leading architecture in the race towards large-
scale quantum computing [1, 2, 3]. However, performance improvements are still
needed in many of the system components, including those used to measure the
states of qubits. To measure a qubit state, an extremely low-power microwave
pulse is used to measure the scattering properties of an electromagnetic
resonator coupled to the qubit [4, 5]. To perform high-fidelity measurements
of the weak microwave signals that carry information about the qubit state, an
exceptionally low-noise amplifier (i.e., close to the standard quantum limit)
with high gain is required [6, 7].

Originally investigated in the late 1960’s [8, 9], interest in Josephson
parametric amplifiers has been renewed due to their potential for meeting
these needs of superconducting qubit technologies [10]. These devices use a
superconducting Josephson junction to facilitate a nonlinear mixing operation
with a strong pump, producing parametric amplification of the desired signal.
To improve the achievable gain, the signal is coupled into the nonlinear
element through a resonator for high interaction time; however, this comes at
the cost of a constrained amplifier bandwidth [7]. Josephson traveling-wave
parametric amplifiers (JTWPAs) were proposed as a wideband alternative, but
their lackluster gain combined with challenging manufacturing requirements
made them difficult to implement [11]. While these devices can achieve long
interaction times by embedding thousands of Josephson junctions in a
transmission line, amplification is inhibited by phase mismatches produced
through pump-induced dispersion [12]. However, modern fabrication techniques
have enabled the development of resonant phase matching (RPM), where resonant
loads are periodically embedded into the transmission line to compensate for
deleterious phase mismatch and obtain significantly higher gain [12, 13]. Over
time, the RPM technique has led JTWPAs to become a promising technology for
enabling qubit readout due to their wide bandwidth, high gain, and ultra low-
noise operation [11, 14, 15, 16]. However, alternative techniques for
achieving phase matching have also been developed, such as periodically
varying the transmission line impedance [17] and modifying the nonlinearity
through flux biasing of Josephson junction-based circuit elements [18].

While modern fabrication technologies and design techniques have allowed high
performance JTWPAs to be realized, their broader application is still limited
by manufacturing variations. Fluctuations between junction properties across a
device can produce impedance mismatches, leading to inconsistent performance
between JTWPAs [19, 20, 21]. Additionally, even small variations in the
resonant loads of RPM architectures can disrupt the delicate phase matching
relationship and significantly reduce gain [22]. Ideally, JTWPA designers
could readily simulate the impact of these variations during the design
process and make adjustments to reduce their severity. However, most JTWPA
models treat the device as a continuum [23, 24, 17, 25] to develop analytical
solutions, preventing variations between elements of the circuit from being
modeled.

To combat this issue, previous approaches have accounted for parameter
variations by introducing a position-dependent wavenumber [22, 26, 12]. In
these methods, analytical expressions are used to describe segments of a JTWPA
with a fixed wavenumber, which are stitched together numerically to simulate
the whole JTWPA. While these approaches can accommodate phase mismatch due to
parameter variations, they do not incorporate impedance mismatches between
segments, which cause reflections that alter performance. As a result, these
methods cannot fully capture the impact that manufacturing variability can
have on a JTWPA design. Aspects of these effects can be captured in numerical
methods using lumped-element circuit simulators such as WRSPICE [27, 28, 29].
However, these approaches, along with the others, are not naturally extensible
to a 3D full-wave method, preventing them from modeling explicit device
geometries. Instead, these methods must estimate circuit parameters from the
explicit geometry, which can often be a user-intensive and potentially
inaccurate procedure that is avoided in a full-wave approach.

In this work, we develop a 1D multiphysics numerical method which treats the
Josephson junctions and transmission lines in a JTWPA as coupled dynamical
systems. The transmission lines are discretized via a finite element time
domain (FETD) approach, while the differential equations describing junction
dynamics are discretized using finite differencing. Interaction between the
dynamical variables is facilitated through a leap-frog time-marching
procedure. In this treatment, individual JTWPA elements can easily be varied
to model arbitrary distributions for the properties of those components.
Additionally, this method inherently considers higher-order nonlinear effects
such as the generation and impact of harmonics. As JTWPA designs evolve,
understanding these effects is increasingly important, but they are difficult
to characterize using analytical approaches [27]. Similarly, pump depletion is
naturally accounted for in this formulation, allowing the dynamic range of a
design to easily be determined. The method can also be extended to a 3D full-
wave description in the future to provide superior accuracy in addition to the
aforementioned benefits.

Preliminary results on this formulation were reported in [30]. This work
significantly expands on [30] by providing a more detailed formulation of the
dynamical equations, as well as extending the method to incorporate RPM.
Further, complete details on the numerical discretization approach are also
presented. Finally, new numerical results are presented, including a case-
study analyzing the impact of process variations in resonant loads and
Josephson junctions.

This work is structured as follows. In Section II, the equations of motion and
boundary conditions for a unit cell of a JTWPA are derived, both with and
without RPM. Following that, Section III develops the 1D hybrid method to
discretize a unit cell, and extends the description to a full JTWPA for each
case. Section IV validates the method using analytical results for each type
of architecture, and demonstrates the effects of variations in the properties
of Josephson junctions and resonant loads on amplification. In Section V, we
discuss the conclusions drawn from these results and the direction for future
development of the method.

##  II Formulation

Figure 1: Schematic of a basic single JTWPA unit cell, consisting of a
Josephson junction with a nonlinear inductance ($L_{\mathrm{J}}$, denoted by
an “X”) in parallel with a linear capacitance ($C_{\mathrm{J}}$).

In this section, the equations of motion (EoMs) and boundary conditions
applicable to a single JTWPA unit cell are derived using techniques from the
Lagrangian mechanics of continuous systems [31]. In Section II-A, this is done
for the “basic JTWPA unit cell” of Fig. 1, while in Section II-B the
derivation is extended to include a resonant load so architectures using RPM
can be analyzed.

###  II-A Basic JTWPA Unit Cell

Deriving EoMs and boundary conditions for the unit cell of Fig. 1 using
Lagrangian mechanics requires a Lagrangian for each portion of the system.
Starting with the transmission lines, the Lagrangian for each line can be
written as

$\displaystyle\mathcal{L}_{\mathrm{TX}_{i}}=\frac{1}{2}\int\big{(}C_{i}V^{2}(z,t)-L_{i}{I}^{2}(z,t)\big{)}\,\mathrm{d}z,$
(1)

where $C_{i}$ and $L_{i}$ are the per-unit-length capacitance and inductance
of the transmission line. However, to derive EoMs the Lagrangian must be
expressed in terms of a generalized coordinate and its derivatives rather than
voltage and current. To accomodate this requirement, we use the node flux
$\phi(z,t)$, which is related to the voltage and current by
$\dot{\phi}(z,t)=V(z,t)$ and $\phi^{\prime}(z,t)=-LI(z,t)$, where
$\dot{\phi}(z,t)\equiv\frac{\partial}{\partial t}\phi(z,t)$ and
$\phi^{\prime}(z,t)\equiv\frac{\partial}{\partial z}\phi(z,t)$ [32]. The
resulting Lagrangian is

$\displaystyle\mathcal{L}_{\mathrm{TX}_{i}}=\frac{1}{2}\int\big{(}C_{i}\dot{\phi}^{2}(z,t)-L_{i}^{-1}{\phi^{\prime}}^{2}(z,t)\big{)}\,\mathrm{d}z.$
(2)

The Josephson junction Lagrangian is developed in [33] as

$\displaystyle\mathcal{L}_{\mathrm{JJ}}=\frac{C_{\mathrm{J}}}{2}\dot{\phi}^{2}_{\mathrm{J}}(t)+\frac{\hbar
I_{\mathrm{c}}}{2e}\cos\bigg{(}\frac{2e}{\hbar}\phi_{\mathrm{J}}(t)\bigg{)},$
(3)

where $e$ is the elementary charge, $\hbar$ is the reduced Planck’s constant,
and $\phi_{\mathrm{J}}(t)\equiv\phi(z_{1},t)-\phi(z_{2},t)$ is the junction
flux. The Josephson junction is characterized by the junction capacitance
$C_{\mathrm{J}}$ and critical current $I_{\mathrm{c}}$, which is related to
the Josephson inductance $L_{\mathrm{J}}$ through
$I_{\mathrm{c}}={\hbar}/(2eL_{\mathrm{J}})$.

To derive the EoMs and boundary conditions in Lagrangian mechanics, we must
first introduce the action functional, $S=\int\mathcal{L}\,\mathrm{d}t$. The
principle of least action states that a physical system will always operate
such that the action is minimized, meaning $\delta S=0$ [31]. Typically, the
principle of least action can be satisfied using the Euler-Lagrange equation.
However, for our continuous systems with non-trivial interconnections as in
Fig. 1, the typical Euler-Lagrange equation is no longer applicable. Hence,
here we instead take the variation of the action functional directly. We
simplify this process by decomposing $\delta S$ into terms for each subsystem,

$\displaystyle\delta S=\delta S_{\mathrm{TX}_{1}}+\delta
S_{\mathrm{TX}_{2}}+\delta S_{\mathrm{JJ}}.$ (4)

For the first transmission line, $\delta S_{\mathrm{TX}_{1}}$ is given by

$\delta
S_{\mathrm{TX}_{1}}=\int\bigg{(}\delta\dot{\phi}(z,t)\frac{\partial\mathcal{L}_{\mathrm{TX}_{1}}}{\partial\dot{\phi}(z,t)}+\delta\phi^{\prime}(z,t)\frac{\partial\mathcal{L}_{\mathrm{TX}_{1}}}{\partial\phi^{\prime}(z,t)}\bigg{)}\,\mathrm{d}t.$
(5)

By evaluating each derivative using (2), (5) becomes

$\delta
S_{\mathrm{TX}_{1}}=\int\int_{0}^{z_{1}}\big{(}C_{1}\dot{\phi}(z,t)\delta\dot{\phi}(z,t)\\\
-L_{1}^{-1}\phi^{\prime}(z,t)\delta\phi^{\prime}(z,t)\big{)}\,\mathrm{d}z\,\mathrm{d}t.$
(6)

To remove derivatives from the variations, we employ integration by parts. For
the temporal variation $\delta\dot{\phi}(z,t)$, it is appropriate to consider
“fixed endpoint” variations so that the boundary terms from the integration by
parts can be neglected [31]. However, for the spatial variation
$\delta\phi^{\prime}(z,t)$, we need to consider variations with “variable
endpoints” to eventually derive a spatial boundary condition at $z=z_{1}$ that
accounts for the effect of the Josephson junction [31]. As a result, (6)
becomes

$\delta
S_{\mathrm{TX}_{1}}=\int\bigg{(}\int_{0}^{z_{1}}\big{(}-C_{1}\ddot{\phi}(z,t)+L_{1}^{-1}\phi^{\prime\prime}(z,t)\big{)}\delta\phi(z,t)\,\mathrm{d}z\\\
-L_{1}^{-1}\phi^{\prime}(z_{1},t)\delta\phi(z_{1},t)\bigg{)}\,\mathrm{d}t,$
(7)

Similarly, for the other transmission line we find $\delta
S_{\mathrm{TX}_{2}}$ is

$\delta
S_{\mathrm{TX}_{2}}=\int\bigg{(}\int_{z_{2}}^{a}\big{(}-C_{2}\ddot{\phi}(z,t)+L_{2}^{-1}\phi^{\prime\prime}(z,t)\big{)}\delta\phi(z,t)\,\mathrm{d}z\\\
+L_{2}^{-1}\phi^{\prime}(z_{2},t)\delta\phi(z_{2},t)\bigg{)}\,\mathrm{d}t.$
(8)

Figure 2: Schematic of a JTWPA unit cell for a RPM design. The added resonant
load and coupling capacitor are outlined in red.

Moving to the Josephson junction, the change of independent variables in (3)
leads to the expression

$\displaystyle\delta
S_{\mathrm{JJ}}=\int\bigg{(}\delta\phi_{\mathrm{J}}(t)\frac{\partial\mathcal{L}}{\partial\phi_{\mathrm{J}}(t)}+\delta\dot{\phi}_{\mathrm{J}}(t)\frac{\partial\mathcal{L}}{\partial\dot{\phi}_{\mathrm{J}}(t)}\bigg{)}\,\mathrm{d}t.$
(9)

By evaluating the derivatives using (3), we obtain

$\delta
S_{\mathrm{JJ}}=\int\bigg{[}C_{\mathrm{J}}\dot{\phi}_{\mathrm{J}}(t)\big{(}\delta\dot{\phi}(z_{1},t)-\delta\dot{\phi}(z_{2},t)\big{)}\\\
-I_{\mathrm{c}}\sin\bigg{(}\frac{2e}{\hbar}\phi_{\mathrm{J}}(t)\bigg{)}\big{(}\delta\phi(z_{1},t)-\delta\phi(z_{2},t)\big{)}\bigg{]}\,\mathrm{d}t,$
(10)

where substitutions were made for $\delta\phi_{\mathrm{J}}(t)$ and
$\delta\dot{\phi}_{\mathrm{J}}(t)$ using
$\phi_{\mathrm{J}}(t)=\phi(z_{1},t)-\phi(z_{2},t)$. Integrating by parts
produces

$\delta
S_{\mathrm{JJ}}=-\int\bigg{[}C_{\mathrm{J}}\ddot{\phi}_{\mathrm{J}}(t)\big{(}\delta\phi(z_{1},t)-\delta\phi(z_{2},t)\big{)}\\\
+I_{\mathrm{c}}\sin\bigg{(}\frac{2e}{\hbar}\phi_{\mathrm{J}}(t)\bigg{)}\big{(}\delta\phi(z_{1},t)-\delta\phi(z_{2},t)\big{)}\bigg{]}\,\mathrm{d}t,$
(11)

which is used alongside (7) and (8) to evaluate (4),

$\delta
S=\int\bigg{[}\int_{0}^{z_{1}}\big{(}-C_{1}\ddot{\phi}(z,t)+L_{1}^{-1}\phi^{\prime\prime}(z,t)\big{)}\delta\phi(z,t)\,\mathrm{d}z\\\
-\bigg{(}L_{1}^{-1}\phi^{\prime}(z_{1},t)+C_{\mathrm{J}}\ddot{\phi}_{\mathrm{J}}(t)+I_{\mathrm{c}}\sin\big{(}\frac{2e}{\hbar}\phi_{\mathrm{J}}(t)\big{)}\bigg{)}\delta\phi(z_{1},t)\\\
+\bigg{(}L_{2}^{-1}\phi^{\prime}(z_{2},t)+C_{\mathrm{J}}\ddot{\phi}_{\mathrm{J}}(t)+I_{\mathrm{c}}\sin\big{(}\frac{2e}{\hbar}\phi_{\mathrm{J}}(t)\big{)}\bigg{)}\delta\phi(z_{2},t)\\\
+\int_{z_{2}}^{a}\big{(}-C_{2}\ddot{\phi}(z,t)+L_{2}^{-1}\phi^{\prime\prime}(z,t)\big{)}\delta\phi(z,t)\,\mathrm{d}z\bigg{]}\mathrm{d}t.$
(12)

The principle of least action states that $\delta S=0$ for all $t$. Since each
variation is nonzero by definition, the coefficient of each variation in (12)
must equate to zero. This step produces the equations of motion and boundary
conditions,

$\displaystyle\phi^{\prime\prime}(z,t)-L_{1}C_{1}\ddot{\phi}(z,t)=0,\quad
0<z<z_{1},$ (13) $\displaystyle
L_{1}^{-1}\phi^{\prime}(z_{1},t)=-C_{\mathrm{J}}\ddot{\phi}_{\mathrm{J}}(t)-I_{\mathrm{c}}\sin\bigg{(}\frac{2e}{\hbar}\phi_{\mathrm{J}}(t)\bigg{)},$
(14) $\displaystyle
L_{2}^{-1}\phi^{\prime}(z_{2},t)=-C_{\mathrm{J}}\ddot{\phi}_{\mathrm{J}}(t)-I_{\mathrm{c}}\sin\bigg{(}\frac{2e}{\hbar}\phi_{\mathrm{J}}(t)\bigg{)},$
(15) $\displaystyle\phi^{\prime\prime}(z,t)-L_{2}C_{2}\ddot{\phi}(z,t)=0,\quad
z_{2}<z<a,$ (16)

Clearly, (13) and (16) describe waves propagating down each transmission line.
The boundary conditions (14) and (15) introduce nonlinearity through
interactions with the Josephson junction. To move towards allowing this
nonlinear set of equations to be solved with a linear time-marching procedure,
we treat $\phi_{\mathrm{J}}(t)$ as a second dynamical variable with EoM

$C_{\mathrm{J}}\ddot{\phi}_{\mathrm{J}}(t)+I_{\mathrm{c}}\sin\bigg{(}\frac{2e}{\hbar}\phi_{\mathrm{J}}(t)\bigg{)}\\\
=-\frac{1}{2L_{1}}\phi^{\prime}(z_{1},t)-\frac{1}{2L_{2}}\phi^{\prime}(z_{2},t).$
(17)

This is obtained by summing (14) and (15). A numerical procedure to enable a
linear solution for time-marching the dynamical variables in (13), (16), and
(17) together will be discussed in Section III.

###  II-B RPM JTWPA Unit Cell

To implement RPM, resonant loads of capacitance $C_{\mathrm{r}}$ and
inductance $L_{\mathrm{r}}$ are periodically coupled into unit cells through
capacitance $C_{\mathrm{c}}$, as shown in Fig. 2. To determine how this change
impacts the equations of motion, we reconstruct the Lagrangian for the second
transmission line, which becomes

$\mathcal{L}_{\mathrm{TX}_{2}}=\frac{1}{2}\int_{z_{2}}^{a}\big{(}C_{2}\dot{\phi}^{2}(z,t)-L_{2}^{-1}{\phi^{\prime}}^{2}(z,t)\big{)}\,\mathrm{d}z\\\
+\frac{1}{2}C_{\mathrm{c}}\big{(}\dot{\phi}(z_{3},t)-\dot{\phi}_{\mathrm{r}}(t)\big{)}^{2}\\!+\frac{1}{2}\big{(}C_{\mathrm{r}}\dot{\phi}_{\mathrm{r}}^{2}(t)-L_{\mathrm{r}}^{-1}{\phi}^{2}_{\mathrm{r}}(t)\big{)},$
(18)

where $\phi_{\mathrm{r}}(t)$ is the flux stored in the resonator. Therefore,
the variation of the action functional is given by

$\delta
S_{\mathrm{TX}_{2}}=\int\bigg{(}\delta\phi_{\mathrm{r}}(z,t)\frac{\partial\mathcal{L}}{\partial\phi_{\mathrm{r}}(z,t)}\\\
+\delta\dot{\phi_{\mathrm{r}}}(z,t)\frac{\partial\mathcal{L}}{\partial\dot{\phi_{\mathrm{r}}}(z,t)}+\delta\phi^{\prime}(z,t)\frac{\partial\mathcal{L}}{\partial\phi^{\prime}(z,t)}\\\
+\delta\dot{\phi}(z,t)\frac{\partial\mathcal{L}}{\partial\dot{\phi}(z,t)}+\delta\dot{\phi}(z_{3},t)\frac{\partial\mathcal{L}}{\partial\dot{\phi}(z_{3},t)}\bigg{)}\,\mathrm{d}t.$
(19)

Substituting (18) into (19) produces the expression

$\delta
S_{\mathrm{TX}_{2}}=\int\bigg{(}\int_{z_{2}}^{a}\big{(}C_{1}\dot{\phi}(z,t)\delta\dot{\phi}(z,t)\\\
-L_{1}^{-1}\phi^{\prime}(z,t)\delta\phi^{\prime}(z,t)\big{)}\,\mathrm{d}z\\\
+C_{\mathrm{r}}\dot{\phi}_{\mathrm{r}}(t)\delta\dot{\phi}_{\mathrm{r}}(t)-L_{\mathrm{r}}^{-1}\phi_{\mathrm{r}}(t)\delta\phi_{\mathrm{r}}(t)\\\
+C_{\mathrm{c}}\big{(}\dot{\phi}(z_{3},t)-\dot{\phi}_{\mathrm{r}}(t)\big{)}\big{(}\delta\dot{\phi}(z_{3},t)-\delta\dot{\phi}_{\mathrm{r}}(t)\big{)}\bigg{)}\,\mathrm{d}t,$
(20)

which after integrating by parts becomes

$\delta
S_{\mathrm{TX}_{2}}=\int\bigg{[}\int_{z_{2}}^{a}\big{(}-C_{2}\ddot{\phi}(z,t)+L_{2}^{-1}\phi^{\prime\prime}(z,t)\big{)}\delta\phi(z,t)\,\mathrm{d}z\\\
+L_{2}^{-1}\phi^{\prime}(z_{2},t)\delta\phi(z_{2},t)\\\
-\bigg{(}C_{\mathrm{r}}\ddot{\phi}_{\mathrm{r}}(t)-L_{\mathrm{r}}^{-1}\phi_{\mathrm{r}}(t)+C_{\mathrm{c}}\big{(}\ddot{\phi}(z_{3},t)-\ddot{\phi}_{\mathrm{r}}(t)\big{)}\bigg{)}\delta\phi_{\mathrm{r}}(t)\\\
-C_{\mathrm{c}}\big{(}\ddot{\phi}(z_{3},t)-\ddot{\phi}_{\mathrm{r}}(t)\big{)}\delta\phi(z_{3},t)\bigg{]}\mathrm{d}t.$
(21)

This expression can be used in conjunction with (7) and (11) to determine
$\delta S$ for the whole RPM unit cell,

$\delta
S=\int\bigg{[}\int_{0}^{z_{1}}\big{(}-C_{1}\ddot{\phi}(z,t)+L_{1}^{-1}\phi^{\prime\prime}(z,t)\big{)}\delta\phi(z,t)\,\mathrm{d}z\\\
+\int_{z_{2}}^{a}\big{(}-C_{2}\ddot{\phi}(z,t)+L_{2}^{-1}\phi^{\prime\prime}(z,t)\big{)}\delta\phi(z,t)\,\mathrm{d}z\\\
-\bigg{(}L_{1}^{-1}\phi^{\prime}(z_{1},t)+\frac{C_{\mathrm{J}}}{2}\ddot{\phi}_{\mathrm{J}}(t)+I_{\mathrm{c}}\sin\big{(}\frac{2e}{\hbar}\phi_{\mathrm{J}}(t)\big{)}\bigg{)}\delta\phi(z_{1},t)\\\
+\bigg{(}L_{2}^{-1}\phi^{\prime}(z_{2},t)+\frac{C_{\mathrm{J}}}{2}\ddot{\phi}_{\mathrm{J}}(t)+I_{\mathrm{c}}\sin\big{(}\frac{2e}{\hbar}\phi_{\mathrm{J}}(t)\big{)}\bigg{)}\delta\phi(z_{2},t)\\\
-\bigg{(}C_{\mathrm{r}}\ddot{\phi}_{\mathrm{r}}(t)-L_{\mathrm{r}}^{-1}\phi_{\mathrm{r}}(t)+C_{\mathrm{c}}\big{(}\ddot{\phi}(z_{3},t)-\ddot{\phi}_{\mathrm{r}}(t)\big{)}\bigg{)}\delta\phi_{\mathrm{r}}(t)\\\
-C_{\mathrm{c}}\big{(}\ddot{\phi}(z_{3},t)-\ddot{\phi}_{\mathrm{r}}(t)\big{)}\delta\phi(z_{3},t)\bigg{]}\mathrm{d}t.$
(22)

Applying the principle of least action to (22) leads to

$\displaystyle\phi^{\prime\prime}(z,t)-L_{1}C_{1}\ddot{\phi}(z,t)=0,\quad
0<z<z_{1},$ (23) $\displaystyle
L_{1}^{-1}\phi^{\prime}(z_{1},t)=-C_{\mathrm{J}}\ddot{\phi}_{\mathrm{J}}(t)-I_{\mathrm{c}}\sin\bigg{(}\frac{2e}{\hbar}\phi_{\mathrm{J}}(t)\bigg{)},$
(24) $\displaystyle
L_{2}^{-1}\phi^{\prime}(z_{2},t)=-C_{\mathrm{J}}\ddot{\phi}_{\mathrm{J}}(t)-I_{\mathrm{c}}\sin\bigg{(}\frac{2e}{\hbar}\phi_{\mathrm{J}}(t)\bigg{)},$
(25) $\displaystyle\phi^{\prime\prime}(z,t)-L_{2}C_{2}\ddot{\phi}(z,t)=0,\quad
z_{2}<z<z_{3},$ (26)
$\displaystyle\phi^{\prime\prime}(z_{3},t)-L_{2}(C_{2}+C_{\mathrm{c}})\ddot{\phi}(z_{3},t)=-L_{2}C_{\mathrm{c}}\ddot{\phi}_{\mathrm{r}}(t),$
(27) $\displaystyle\phi^{\prime\prime}(z,t)-L_{2}C_{2}\ddot{\phi}(z,t)=0,\quad
z_{3}<z<a,$ (28)
$\displaystyle\phi_{\mathrm{r}}(t)-L_{\mathrm{r}}(C_{\mathrm{r}}+C_{\mathrm{c}})\ddot{\phi}_{\mathrm{r}}(t)=-L_{\mathrm{r}}C_{\mathrm{c}}\ddot{\phi}(z_{3},t).$
(29)

Since (24) and (25) are unchanged from (14) and (15), (17) still holds for the
RPM unit cell. A numerical procedure to enable a linear solution for time-
marching the dynamical variables in (17), (23), and (26) to (29) together will
be discussed in Section III.

##  III Discretization

Here, we detail the development of the hybrid model from the EoMs for each
unit cell. In this method, 1D FETD is applied to discretize $\phi(z,t)$ along
the transmission lines. FETD is chosen to simplify the transition to a 3D
full-wave implementation in the future. Temporal discretization of
$\phi_{\mathrm{J}}(t)$ is performed using central differencing, which will be
demonstrated to produce a linear time-marching equation. The direct
relationship between the dynamical variables is eliminated using leap-frog
time-marching, which allows $\phi(z,t)$ to be evaluated without the need for a
nonlinear solver.

This section is structured as follows. First, the method for the basic JTWPA
unit cell is developed in Section III-A by deriving time-marching equations
for the dynamical variables. Additionally, we introduce the leap-frog time-
stepping procedure to linearly advance the dynamical variables in time. In
Section III-B, the unit cells are combined and integrated with driving and
terminating components to model a full JTWPA. Finally, the time-marching
equations are adjusted to accomodate the addition of resonant loads to model a
RPM JTWPA in Section III-C.

###  III-A Basic JTWPA Unit Cell

To begin, we discretize the node flux $\phi(z,t)$ along the transmission
lines, which is described by (13) to (16). We follow a standard FETD process
outlined in [34, 32] with first-order (triangular) nodal basis functions
$N_{i}(z)$ such that

$\displaystyle\phi(z,t)=\sum_{i}\phi_{i}(t)N_{i}(z).$ (30)

This process leads to a matrix equation of the form

$\displaystyle[T_{\mathrm{u}}]\frac{d^{2}}{dt^{2}}\\{\phi(t)\\}+[S_{\mathrm{u}}]\\{\phi(t)\\}=\\{f_{\mathrm{u}}(t)\\},$
(31)

where $\\{\phi(t)\\}_{i}=\phi_{i}(t)$ and

$[T_{\mathrm{u}}]_{ij}=L_{1}C_{1}\int_{0}^{z_{1}}N_{i}(z)N_{j}(z)\,\mathrm{d}z\\\
+L_{1}C_{\mathrm{J}}\big{(}\delta_{jM}-\delta_{j(M+1)}\big{)}\delta_{iM}+L_{2}C_{2}\int_{z_{2}}^{a}N_{i}(z)N_{j}(z)\,\mathrm{d}z\\\
-L_{2}C_{\mathrm{J}}\big{(}\delta_{jM}-\delta_{j(M+1)}\big{)}\delta_{i(M+1)},$
(32)
$\displaystyle[S_{\mathrm{u}}]_{ij}=\int_{0}^{z_{1}}N_{i}^{\prime}(z)N_{j}^{\prime}(z)\,\mathrm{d}z+\int_{z_{2}}^{a}N_{i}^{\prime}(z)N_{j}^{\prime}(z)\,\mathrm{d}z,$
(33)
$\displaystyle\\{f_{\mathrm{u}}(t)\\}_{i}=I_{\mathrm{c}}\big{(}L_{2}\delta_{i(M+1)}-L_{1}\delta_{iM}\big{)}\sin\bigg{(}\frac{2e}{\hbar}\phi_{\mathrm{J}}(t)\bigg{)},$
(34)

where $\delta_{ij}$ is a Kronecker delta, and nodes $M$ and $M+1$ are
positioned at $z_{1}$ and $z_{2}$ respectively. For the junction flux,
applying (30) leads (17) to become

$C_{\mathrm{J}}\ddot{\phi}_{\mathrm{J}}(t)+I_{\mathrm{c}}\sin\bigg{(}\frac{2e}{\hbar}\phi_{\mathrm{J}}(t)\bigg{)}\\\
=\sum_{i}\bigg{(}-\frac{1}{2L_{1}}N_{i}^{\prime}(z_{1})-\frac{1}{2L_{2}}N_{i}^{\prime}(z_{2})\bigg{)}\phi_{i}(t)\\\
=-\frac{1}{2L_{1}\Delta z}\big{(}\phi_{M}(t)-\phi_{M-1}(t)\big{)}\\\
-\frac{1}{2L_{2}\Delta z}\big{(}\phi_{M+2}(t)-\phi_{M+1}(t)\big{)},$ (35)

where the final equality was reached by taking piecewise derivatives of the
triangular basis functions.

At this point, $\\{\phi(t)\\}$ and $\phi_{\mathrm{J}}(t)$ are still directly
related through $\phi_{\mathrm{J}}(t)\equiv\phi_{M}(t)-\phi_{M+1}(t)$, which
invalidates the source treatment of $\phi_{\mathrm{J}}(t)$ in (34). However,
leap-frogging between the dynamical variables in the time discretization can
be exploited to circumvent this issue. More explicitly, by solving for each
dynamical variable at alternating half-timesteps, the direct relationship
between them is no longer enforced, enabling the “source treatment” for
$\phi_{\mathrm{J}}(t)$ in (34). Applying this method using central
differencing yields time-marching equations for each dynamical variable,

$[T_{\mathrm{u}}]\\{\phi^{(n+1)}\\}=\frac{(\Delta
t)^{2}}{2}\big{(}\\{f_{\mathrm{u}}^{(n+1/2)}\\}+\\{f_{\mathrm{u}}^{(n-1/2)}\\}\big{)}\\\
-\big{(}(\Delta
t)^{2}[S_{\mathrm{u}}]-2[T_{\mathrm{u}}]\big{)}\\{\phi^{(n)}\\}-[T_{\mathrm{u}}]\\{\phi^{(n-1)}\\}$
(36)

and

$\phi_{\mathrm{J}}^{(n+3/2)}=2\phi_{\mathrm{J}}^{(n+1/2)}-\frac{I_{\mathrm{c}}(\Delta
t)^{2}}{C_{\mathrm{J}}}\sin\bigg{(}\frac{2e}{\hbar}\phi_{\mathrm{J}}^{(n+1/2)}\bigg{)}\\\
-\phi_{\mathrm{J}}^{(n-1/2)}-\frac{(\Delta t)^{2}}{4L_{1}C_{\mathrm{J}}\Delta
z}\big{(}\phi_{M}^{(n+1)}+\phi_{M}^{(n)}-\phi_{M-1}^{(n+1)}-\phi_{M-1}^{(n)}\big{)}\\\
-\frac{(\Delta t)^{2}}{4L_{2}C_{\mathrm{J}}\Delta
z}\big{(}\phi_{M+2}^{(n+1)}+\phi_{M+2}^{n}-\phi_{M+1}^{(n+1)}-\phi_{M+1}^{n}\big{)},$
(37)

where $\phi_{M}^{(n)}$ is the flux at node $M$ and timestep $n$. In
formulating (36) and (37), averaging was used to approximate unknown values,
i.e.:
$\phi_{M}^{(n+1/2)}\approx\frac{1}{2}\big{(}\phi_{M}^{(n+1)}+\phi_{M}^{(n)}\big{)}$.
Notably, using central differencing to discretize (35) led to a time-marching
equation for $\phi_{\mathrm{J}}(t)$ that can be solved linearly, despite the
nonlinearity of (35). As a result, both dynamical variables can be evolved in
time without the need for a nonlinear solver.

###  III-B Full JTWPA Simulation

Figure 3: Schematic of a JTWPA. Capacitive loads are introduced at the source
and termination of the JTWPA to aid in filtering out low-frequency spurious
modes that lead to late-time instability.

The method is extended from a single unit cell to a JTWPA by connecting $N$
unit cells with per-unit-length inductance $L_{\mathrm{u}_{n}}$ and
capacitance $C_{\mathrm{u}_{n}}$ together at their endpoints, as shown in Fig.
3. To drive the JTWPA, a voltage source $V_{\mathrm{s}}(t)$ with resistance
$R_{\mathrm{s}}$ is connected to the the input side through a transmission
line with per-unit-length inductance $L_{\mathrm{in}}$ and capacitance
$C_{\mathrm{in}}$. On the output side, a load $R_{\mathrm{l}}$ is connected
through another transmission line characterized by per-unit-length parameters
$L_{\mathrm{out}}$ and $C_{\mathrm{out}}$ to terminate the simulation region.
The source and load capacitors $C_{\mathrm{s}}$ and $C_{\mathrm{l}}$ are
introduced to aid in filtering out low-frequency spurious modes that lead to
late-time instability, which will be discussed in more detail in Section IV.
The boundary conditions for the source and load are introduced in [32] as

$\displaystyle\phi^{\prime}(z_{\mathrm{s}},t)=-\frac{L_{\mathrm{in}}}{R_{\mathrm{s}}}[V_{\mathrm{s}}(t)-\dot{\phi}(z_{\mathrm{s}},t)],$
(38)
$\displaystyle\phi^{\prime}(z_{\mathrm{l}},t)=-\frac{L_{\mathrm{out}}}{R_{\mathrm{l}}}\dot{\phi}(z_{\mathrm{l}},t),$
(39)

where $z_{\mathrm{s}}$ and $z_{\mathrm{l}}$ are the positions of the source
and load. Since these terms contain a first derivative in time, the matrix
equation of (31) needs to be modified to the form

$\displaystyle[T]\frac{d^{2}}{dt^{2}}\\{\phi(t)\\}+[R]\frac{d}{dt}\\{\phi(t)\\}+[S]\\{\phi(t)\\}=\\{f(t)\\}.$
(40)

The matrices for the full JTWPA in this expression become

$\displaystyle[T]=[T_{\mathrm{in}}]+\sum_{n}[T_{\mathrm{u}_{n}}]+[T_{\mathrm{out}}],$
(41)
$\displaystyle[R]_{ij}=\frac{L_{\mathrm{in}}}{R_{\mathrm{s}}}\delta_{i1}\delta_{j1}+\frac{L_{\mathrm{out}}}{R_{l}}\delta_{iN_{e}}\delta_{iN_{e}},$
(42)
$\displaystyle[S]=[S_{\mathrm{in}}]+\sum_{n}[S_{\mathrm{u}_{n}}]+[S_{\mathrm{out}}],$
(43)
$\displaystyle\\{f(t)\\}_{i}=\frac{L_{\mathrm{in}}}{R_{\mathrm{s}}}V_{\mathrm{s}}(t)\delta_{i1}+\sum_{n}\\{f_{\mathrm{u}_{n}}(t)\\},$
(44)

where node $1$ occurs at the source, and node $N_{e}$ is positioned at the
load. Further, $[T_{\mathrm{u}_{n}}]$, $[S_{\mathrm{u}_{n}}]$, and
$\\{f_{\mathrm{u}_{n}}(t)\\}$ refer to (32), (33), and (34) for unit cell $n$.
The input matrices have the form

$[T_{\mathrm{in}}]_{ij}=L_{\mathrm{in}}C_{\mathrm{in}}\int_{z_{\mathrm{s}}}^{z_{0}}N_{i}(z)N_{j}(z)\,\mathrm{d}z\\\
+L_{\mathrm{in}}C_{\mathrm{s}}\big{(}\delta_{11}-\delta_{12}-\delta_{21}+\delta_{22}\big{)},$
(45)
$\displaystyle[S_{\mathrm{in}}]_{ij}=\int_{z_{\mathrm{s}}}^{z_{0}}N_{i}^{\prime}(z)N_{j}^{\prime}(z)\,\mathrm{d}z,$
(46)

where $z_{0}$ is the starting position of the first unit cell. The output
matrices can be expressed as

$[T_{\mathrm{out}}]_{ij}=L_{\mathrm{out}}C_{\mathrm{out}}\int_{z_{\mathrm{e}}}^{z_{\mathrm{l}}}N_{i}(z)N_{j}(z)\,\mathrm{d}z\\\
+L_{\mathrm{out}}C_{\mathrm{l}}\big{(}\delta_{(N_{e}-1)(N_{e}-1)}-\delta_{(N_{e}-1)N_{e}}\\\
-\delta_{N_{e}(N_{e}-1)}+\delta_{N_{e}N_{e}}\big{)},$ (47)
$\displaystyle[S_{\mathrm{out}}]_{ij}=\int_{z_{\mathrm{e}}}^{z_{\mathrm{l}}}N_{i}^{\prime}(z)N_{j}^{\prime}(z)\,\mathrm{d}z,$
(48)

where $z_{\mathrm{e}}$ is the end position of the final unit cell.

With the matrices in (40) constructed, central differencing can be applied to
produce a time-marching equation,

$\bigg{(}[T]+\frac{\Delta t}{2}[R]\bigg{)}\\{\phi^{(n+1)}\\}=(\Delta
t)^{2}\\{f^{(n)}\\}\\\ +\frac{(\Delta
t)^{2}}{2}\big{(}\\{f^{(n+1/2)}\\}+\\{f^{(n-1/2)}\\}\big{)}\\\ -\big{(}(\Delta
t)^{2}[S]-2[T]\big{)}\\{\phi^{(n)}\\}-\bigg{(}[T]-\frac{\Delta
t}{2}[R]\bigg{)}\\{\phi^{(n-1)}\\}.$ (49)

The time marching equation for $\phi_{\mathrm{J}}(t)$ in (37) is unchanged for
the full JTWPA, but now must be independently solved for each Josephson
junction in the circuit.

###  III-C RPM JTWPA

The inclusion of resonant loads introduces the resonator flux
$\phi_{\mathrm{r}}(t)$ as a third dynamical variable for an RPM JTWPA.
However, since (29) is linear, there is no need to incorporate leap-frogging
to evolve $\phi_{\mathrm{r}}(t)$ in time. Instead, the resonator flux
$\phi_{\mathrm{r}}(t)$ can be solved simultaneously with $\phi(z,t)$. To
accomodate this change, we modify (40) to the form

$\begin{bmatrix}T_{\ell\ell}&T_{\ell\mathrm{r}}\\\
T_{\mathrm{r}\ell}&T_{\mathrm{r}\mathrm{r}}\end{bmatrix}\frac{d^{2}}{dt^{2}}\begin{Bmatrix}\phi(t)\\\
\phi_{\mathrm{r}}(t)\end{Bmatrix}+\begin{bmatrix}R_{\ell\ell}&0\\\ 0&0\\\
\end{bmatrix}\frac{d}{dt}\begin{Bmatrix}\phi(t)\\\
\phi_{\mathrm{r}}(t)\end{Bmatrix}\\\ +\begin{bmatrix}S_{\ell\ell}&0\\\
0&S_{\mathrm{r}\mathrm{r}}\end{bmatrix}\begin{Bmatrix}\phi(t)\\\
\phi_{\mathrm{r}}(t)\end{Bmatrix}=\begin{Bmatrix}f(t)\\\ 0\end{Bmatrix},$ (50)

where $\\{\phi_{\mathrm{r}}(t)\\}_{i}=\phi_{\mathrm{r},i}(t)$, or the flux in
the resonant load of the $i$th unit cell. Further, $[R]_{\ell\ell}$,
$[S]_{\ell\ell}$, and $\\{f(t)\\}$ are equivalent to their counterparts in
(42) to (44) for the full JTWPA without RPM. From (23) to (29), the remaining
matrices can be expressed as

$\displaystyle[T_{\ell\ell}]_{ij}=[T]_{ij}+\sum_{n}L_{\mathrm{u}_{n}}C_{\mathrm{c}_{n}}\delta_{iR_{n}}\delta_{jR_{n}}$
(51)
$\displaystyle[T_{\ell\mathrm{r}}]_{ij}=\sum_{n}-L_{\mathrm{u}_{n}}C_{\mathrm{c}_{n}}\delta_{iR_{n}}\delta_{jn},$
(52)
$\displaystyle[T_{\mathrm{r}\ell}]_{ij}=\sum_{n}-L_{\mathrm{r}_{n}}C_{\mathrm{c}_{n}}\delta_{in}\delta_{jR_{n}},$
(53)
$\displaystyle[T_{\mathrm{r}\mathrm{r}}]_{ij}=\sum_{n}L_{\mathrm{r}_{n}}(C_{\mathrm{c}_{n}}+C_{\mathrm{r}_{n}})\delta_{in}\delta_{jn},$
(54) $\displaystyle[S_{\mathrm{r}\mathrm{r}}]=\sum_{n}\delta_{in}\delta_{jn},$
(55)

where each node $R_{n}$ is where the resonant load is located in the $n$th
unit cell, and $L_{\mathrm{r}_{n}}$, $C_{\mathrm{c}_{n}}$, and
$C_{\mathrm{r}_{n}}$ are $L_{\mathrm{r}}$, $C_{\mathrm{c}}$, and
$C_{\mathrm{r}}$ for that unit cell.

The time-marching expression for the RPM JTWPA is the same as the equation for
the basic JTWPA in (49), but with the matrices and vectors modified to the
block form of (50). The time-marching equation for each $\phi_{\mathrm{J}}(t)$
in (37) is unchanged by the addition of resonant loads.

##  IV Results

Here, we present results to validate the numerical method and exhibit its
utility. In Section IV-A, the gain of a basic JTWPA design is calculated from
the results of the numerical method and compared to analytical results. This
procedure is extended to a JTWPA architecture using RPM in Section IV-B, which
also includes a discussion of the strategy used to inhibit late-time
instability. In Section IV-C, the method is applied to simulate the effect of
variations in Josephson junction area across a JTWPA. Further, the resonant
loads are varied in Section IV-D to demonstrate another potential application.

###  IV-A Basic JTWPA Gain Curve

Figure 4: Gain of the architecture of Table I calculated using the analytical
method of [25] and the numerical method of this work. Points between
$5.3\,\mathrm{GHz}$ and $6.6\,\mathrm{GHz}$ are interpolated using a quadratic
fit.

TABLE I: Parameters for the basic JTWPA architecture from [12] Parameter |  Value
---|---
Junction critical current ($I_{\mathrm{c}}$) |  $3.29\,\upmu\mathrm{A}$
Josephson capacitance ($C_{\mathrm{J}}$) |  $329\,\mathrm{fF}$
Per-unit-length capacitance ($C_{\mathrm{u}}$) |  $3.9\,\mathrm{nF/m}$
Unit cell length ($a$) |  $10\,\upmu\mathrm{m}$
Number of unit cells ($N$) |  $2000$
Pump current magnitude ($I_{\mathrm{p}}$) |  $1.645\,\upmu\mathrm{A}$
Pump frequency ($f_{\mathrm{p}}$) |  $5.970\,\mathrm{GHz}$

The method is initially validated by comparing to the analytical results from
[25] for the basic JTWPA design presented in [12]. The specifications for this
JTWPA are detailed in Table I. Notably, the inductance of each unit cell in
this design is assumed to be dominated by $L_{\mathrm{J}}$, meaning
$L_{\mathrm{u}}\approx 0$. This assumption is not suitable for our numerical
method, so a small per-unit-length inductance of
$L_{\mathrm{u}}=1\,\upmu\mathrm{H/m}$ is added to the transmission lines. This
small $L_{\mathrm{u}}$ was found to have a minimal impact on the numerical
results, while also not modifying the FETD stability constraint to a
problematic extent. To match the input and output transmission lines to the
JTWPA, the contributions of the Josephson junctions need to be accounted for,
so $L_{\mathrm{in}}=L_{\mathrm{out}}=L_{\mathrm{u}}+L_{\mathrm{J}}/a$.

The JTWPA is driven with a voltage source

$\displaystyle
V_{\mathrm{s}}(t)=2I_{\mathrm{p}}Z_{0}W(t)\big{(}\sin(\omega_{\mathrm{p}}t)+\alpha\sin(\omega_{\mathrm{s}}t)\big{)},$
(56)

where $\omega_{\mathrm{s}}$ is the signal frequency, $\omega_{\mathrm{p}}$ is
the pump frequency, and
$Z_{0}=(aL_{\mathrm{u}}+L_{\mathrm{J}})/aC_{\mathrm{u}}$ is the intrinsic
impedance of the transmission line. The ratio between the signal and pump
amplitudes is determined by $\alpha$, which is set to $\alpha=1/20$ to prevent
pump depletion from compressing gain. Further, $W(t)$ is an exponential
tapering function with a total width of $50\,\mathrm{ns}$. While analytical
models typically assume the pump and signal are monochromatic, due to our
method operating in the time domain a finite bandwidth is necessary to keep
the total simulation duration at a reasonable level.

The gain $G$ is calculated from simulation results via

$\displaystyle
G=10\log\bigg{(}\frac{\int_{f_{\mathrm{L}}}^{f_{\mathrm{H}}}|\phi_{\mathrm{o}}(f)|^{2}\,\mathrm{d}f}{\int_{f_{\mathrm{L}}}^{f_{\mathrm{H}}}|\phi_{\mathrm{i}}(f)|^{2}\,\mathrm{d}f}\bigg{)},$
(57)

where $\phi_{\mathrm{i}}(f)$ and $\phi_{\mathrm{o}}(f)$ are the Fourier
transform of the node flux at the input and output respectively. The
integration bandwidth is determined by $f_{\mathrm{L}}$ and $f_{\mathrm{H}}$,
which are set to $f_{\mathrm{L}}=f_{\mathrm{s}}-0.5\,\mathrm{GHz}$ and
$f_{\mathrm{H}}=f_{\mathrm{s}}+0.5\,\mathrm{GHz}$. Wideband integration is
necessary due to self-mixing of the signal, which spreads energy from its
center frequency over a wider frequency range. In an experimental setting, the
signal is typically several orders of magnitude smaller than the critical
current, which limits this self-mixing process. In the numerical method, this
process becomes relevant due to the relatively large signal-to-pump ratio,
which was chosen to allow an increased pump bandwidth that helps reduce the
total simulation duration to keep the method more efficient.

Using (57), the gain is extracted from numerical simulations for the
architecture of Table I and compared to the analytical result of [25] in Fig.
4. At signal frequencies between $5.3\,\mathrm{GHz}$ and $6.6\,\mathrm{GHz}$,
the integration bandwidth overlaps with the pump and induces error in the gain
calculation. To allow the gain to be extracted in this region, a quadratic fit
is used to interpolate it from surrounding points. While interpolation is
necessary for this topology, alternative JTWPA designs using three-wave mixing
remove the pump from the operating bandwidth, allowing this process to be
avoided. These architectures will be discussed in more detail in Section V and
considered in future work.

###  IV-B RPM JTWPA Gain Curve

TABLE II: Parameters for the RPM JTWPA architecture from [12] Parameter |  Value
---|---
Junction critical current ($I_{\mathrm{c}}$) |  $3.29\,\upmu\mathrm{A}$
Josephson capacitance ($C_{\mathrm{J}}$) |  $329\,\mathrm{fF}$
Per-unit-length capacitance ($C_{\mathrm{u}}$) |  $3.9\,\mathrm{nF/m}$
Unit cell length ($a$) |  $10\,\upmu\mathrm{m}$
Resonant load coupling capacitance ($C_{\mathrm{c}}$) |  $10\,\mathrm{fF}$
Resonant load capacitance ($C_{\mathrm{r}}$) |  $7.036\,\mathrm{pF}$
Resonant load inductance ($L_{\mathrm{r}}$) |  $100\,\mathrm{pH}$
Number of unit cells ($N$) |  $2000$
Pump current magnitude ($I_{\mathrm{p}}$) |  $1.645\,\upmu\mathrm{A}$
Pump frequency ($f_{\mathrm{p}}$) |  $5.970\,\mathrm{GHz}$
Figure 5: Gain of the architecture of Table II calculated using the analytical
method of [25] and the numerical method of this work. Points between
$5.3\,\mathrm{GHz}$ and $6.6\,\mathrm{GHz}$ are interpolated using a fourth-
order polynomial fit.

The method is also validated using the RPM JTWPA design in [12], for which the
relevant parameter values are listed in Table II. For this architecture, some
additional changes to the voltage source described in (56) must be made. First
of all, the addition of the resonant loads shifts the average intrinsic
impedance of the JTWPA to
$Z_{0}=(aL_{\mathrm{u}}+L_{\mathrm{J}})/(aC_{\mathrm{u}}+C_{\mathrm{c}})$. To
compensate for this, the per-unit-length capacitance of the input and output
transmission lines must also be adjusted to
$C_{\mathrm{in}}=C_{\mathrm{out}}=C_{\mathrm{u}}+C_{\mathrm{c}}/a$.
Additionally, the increased gain of this design requires that the signal-to-
pump ratio is reduced to $\alpha=1/50$ to minimize the effect of pump
depletion. Finally, the pump bandwidth must be reduced for this architecture
due to the addition of the resonant loads. To correct for mismatch between the
signal and pump phases, the resonant loads introduce a pole at
$6\,\mathrm{GHz}$, where its proximity to the pump leads to a phase shift.
However, if the pump is too wideband, this proximity means the sidebands of
the pump will overlap with the pole, resulting in deformation. In an
experimental setting, the pump is typically narrowband enough that this is not
an issue. However, the shortened envelope used in the numerical method
increases the pump bandwidth, causing it to extend into the pole. To minimize
this issue, the width of the tapering function $W(t)$ was increased to
$250\,\mathrm{ns}$.

While the narrowband pump is necessary in this architecture, increasing the
duration of $W(t)$ requires an equivalent increase in total simulation time.
As a result, low-frequency spurious solutions are given time to grow and
eventually produce instability. To address this issue, we applied the
correction method detailed in [35]. This method periodically isolates and
removes the electrostatic portion of the solution to prevent buildup of
spurious solutions. To separate this electrostatic component, a basis for the
nullspace of $[S]$ must be known. For the 1D system, each vector in the basis
is constant over a single segment of transmission line terminated by Josephson
junctions or capacitors, and zero elsewhere.

(a)

(b)

(c)

(d)

Figure 6: Gain of the Table I JTWPA with variations in Josephson junction
area. Each curve is shaded between the minimum and maximum outcomes of 50
simulations. The area of each junction in the JTWPA is independently selected
on a normal distribution with standard deviation $\sigma_{A}$. The mean of the
distribution is (a) constant, or linearly varied by (b) $0.05A$, (c) $0.10A$,
or (d) $0.15A$ across the JTWPA.

Importantly, the correction method of [35] is designed to be performed on
lossless systems. While the JTWPA itself is lossless, $R_{\mathrm{s}}$ and
$R_{\mathrm{l}}$ introduce loss which can produce inaccuracy in the correction
process. To avoid needing to apply the method in these regions, capacitors
$C_{\mathrm{s}}$ and $C_{\mathrm{l}}$ are introduced. These capacitors
separate the source and load from the input and output transmission lines,
which removes them from the corresponding nullspace vectors. As a result, the
correction method can be applied to the input and output transmission lines
without also applying it to the source and load. To prevent instability from
building up in the uncorrected regions and spreading to the JTWPA, the
capacitors are sized at $C_{\mathrm{s}}=C_{\mathrm{l}}=100\,\mathrm{pF}$. At
this value, they present a high impedance to spurious solutions to inhibit
their growth. However, the impedance at the operating frequency range is small
enough that the capacitors do not have a significant impact on the accuracy of
the method.

The gain extrapolated from numerical results for the Table II JTWPA using (57)
is compared against the analytical result from [25] in Fig. 4. For this
architecture, a fourth-order polynomial interpolating function was found to
provide good agreement with the analytical results in the vicinity of the
pump. At $6\,\mathrm{GHz}$, the pole introduced by the resonant loads produces
a large decrease in gain. Due to the use of interpolation in this region, this
effect is not modeled by the numerical method.

###  IV-C Josephson Junction Area Variations

To demonstrate a potential application of the numerical method, we
investigated the effect of variations in Josephson junction area on JTWPA
gain. The junction area $A$ is related to the Josephson inductance and
capacitance by $L_{\mathrm{J}}\propto A$ and $C_{\mathrm{J}}\propto 1/A$ [36].
Fig. 6LABEL:sub@subfig:noRPM-areaGrad=0 demonstrates the result of
independently varying the area of each Josephson junction in the Table I JTWPA
on a normal distribution with mean $A$ and standard deviation $\sigma_{A}$.
For $\sigma_{A}=0.03A$, the impact is relatively minor, and gain remains close
to the analytical result. The differences are more pronounced for
$\sigma_{A}=0.06A$ and $\sigma_{A}=0.09A$, especially at signal frequencies
greater than the pump frequency, where substantial gain ripples occur.

In addition to entirely random variations, significant gradients in Josephson
junction area across a JTWPA have been observed during device manufacturing
[37]. To model these gradients, we varied the mean of the normal distribution
as a linear function of position along the JTWPA. In Fig.
6LABEL:sub@subfig:noRPM-areaGrad=5, the gradient is set so that the mean area
of the final Josephson junction in the JTWPA is $0.05A$ higher than the first
junction. In Fig. 6LABEL:sub@subfig:noRPM-areaGrad=10 and Fig.
6LABEL:sub@subfig:noRPM-areaGrad=15, this procedure was repeated with total
gradients of $0.10A$ and $0.15A$ respectively. These results indicate that
increasing the gradient reduces gain across the bandwidth of the device.

Figure 7: Gain of the Table II JTWPA with variations in Josephson junction
area. Each curve is shaded between the minimum and maximum outcomes of 30
simulations. The area of the Josephson junctions are selected from a normal
distribution with standard deviation $\sigma_{A}$, and with its mean linearly
varied by $0.05A$ across the JTWPA.

The impact of junction area fluctuations was also tested for the JTWPA with
RPM in Fig. 7, where a gradient of $0.05A$ was used. In this case, the
presence of the resonant loads leads to additional gain ripples near the pump
frequency.

###  IV-D Resonant Load Variations

In addition to the Josephson junction area, the method was applied to
investigate how small variations in the values of the resonant load components
can affect JTWPA gain. While these components do not suffer from the same
manufacturing inconsistencies as Josephson junctions, the delicate nature of
the phase matching condition means small changes can severely impact
performance. To demonstrate this, we varied $L_{\mathrm{r}}$ and
$C_{\mathrm{r}}$ for each load on normal distributions with standard
deviations $\sigma_{L_{\mathrm{r}}}$ and $\sigma_{C_{\mathrm{r}}}$. The mean
value for each distribution matches the component value used in Table II. The
result of applying these variations are displayed in Fig. 8 for several values
of $\sigma_{L_{\mathrm{r}}}$ and $\sigma_{C_{\mathrm{r}}}$. For
$\sigma_{L_{\mathrm{r}}(C_{\mathrm{r}})}=0.0010L_{\mathrm{r}}(C_{\mathrm{r}})$
and
$\sigma_{L_{\mathrm{r}}(C_{\mathrm{r}})}=0.0015L_{\mathrm{r}}(C_{\mathrm{r}})$,
the effect of these variations is mostly negligible, indicating that changes
in the pump phase due to the shifting resonance roughly cancel across the
JTWPA. However, the impact at
$\sigma_{L_{\mathrm{r}}(C_{\mathrm{r}})}=0.0020L_{\mathrm{r}}(C_{\mathrm{r}})$,
is clearly much more significant. At this value, the standard deviation is
large enough that the pole may be shifted so that it overlaps with the pump
bandwidth. As a result, the pump is partially absorbed by the resonant load,
leading to the gain drop-off observed in the plot.

##  V Conclusion

Figure 8: Gain of the Table II JTWPA with variations in the resonant loads.
Each curve is shaded between the minimum and maximum outcomes of 30
simulations. For each load, $L_{\mathrm{r}}$ and $C_{\mathrm{r}}$ are
independently chosen from normal distributions with standard deviations
$\sigma_{L_{\mathrm{r}}}$ and $\sigma_{C_{\mathrm{r}}}$, and mean values
matching those in Table II.

In this work, a numerical method for modeling JTPWAs was developed to overcome
the limitations of analytical models. Rather than using a continuum
approximation as typically done in analytical approaches, the Josephson
junctions and transmission lines were treated as coupled subsystems. As a
result, components can easily be individually varied to examine the impact of
manufacturing tolerances on device performance. To linearize the time-marching
equations for the dynamical variables, we employed a leap-frog time-marching
procedure. Gain was extracted from simulations of JTWPA designs with and
without RPM and compared to analytical results to validate the method.
Further, the impact of variations in Josephson junctions and resonant load
components was investigated using the method.

Future work will expand the method to support additional JTWPA architectures
and explore other effects not captured by analytical models. Modern JTWPA
designs integrate Josephson junction-based circuit elements such as
superconducting quantum interference devices (SQUIDs) [17] and superconducting
nonlinear asymmetric inductive elements (SNAILs) [18] to leverage their
additional degrees of freedom. These components can also be used to develop
three-wave mixing JTWPAs which are inherently phase matched [24], and can use
a flux-driven topology to achieve high dynamic range [37, 38]. However, the
performance of these JTWPAs is heavily impacted by the generation of
harmonics, an effect which is not captured by analytical models [27]. By
incorporating components such as SQUIDs and SNAILs into our multiphysics
method, this issue can be explored in greater detail.

In addition to supporting more JTWPA architectures, future work will extend
the numerical method to a 3D full-wave description. While analytical device
models are inherently limited to 1D, the need to extract circuit parameters of
JTWPA geometries limits the accuracy of those methods and adds a significant
pre-processing step. Meanwhile, a 3D numerical method would be capable of
modeling the intended and non-ideal behavior of JTWPAs with unprecedented
accuracy.

## References

  * [1] F. Arute, K. Arya, R. Babbush, D. Bacon, J. Bardin, R. Barends, R. Biswas, S. Boixo, F. Brandao, D. Buell _et al._ , “Quantum supremacy using a programmable superconducting processor,” _Nature_ , vol. 574, p. 505–510, Oct 2019.
  * [2] Y. Wu, W.-S. Bao, S. Cao, F. Chen, M.-C. Chen, X. Chen, T.-H. Chung, H. Deng, Y. Du, D. Fan _et al._ , “Strong quantum computational advantage using a superconducting quantum processor,” _Physical Review Letters_ , vol. 127, no. 18, p. 180501, 2021.
  * [3] P. Jurcevic, A. Javadi-Abhari, L. S. Bishop, I. Lauer, D. F. Bogorin, M. Brink, L. Capelluto, O. Günlük, T. Itoko, N. Kanazawa _et al._ , “Demonstration of quantum volume 64 on a superconducting quantum computing system,” _Quantum Science and Technology_ , vol. 6, no. 2, p. 025020, 2021\.
  * [4] T. E. Roth, R. Ma, and W. C. Chew, “The transmon qubit for electromagnetics engineers: An introduction,” _IEEE Antennas and Propagation Magazine_ , vol. 65, no. 2, pp. 8–20, 2023.
  * [5] P. Krantz, M. Kjaergaard, F. Yan, T. P. Orlando, S. Gustavsson, and W. D. Oliver, “A quantum engineer’s guide to superconducting qubits,” _Applied physics reviews_ , vol. 6, no. 2, 2019.
  * [6] T. Walter, P. Kurpiers, S. Gasparinetti, P. Magnard, A. Potočnik, Y. Salathé, M. Pechal, M. Mondal, M. Oppliger, C. Eichler, and A. Wallraff, “Rapid high-fidelity single-shot dispersive readout of superconducting qubits,” _Phys. Rev. Appl._ , vol. 7, p. 054020, May 2017.
  * [7] J. Aumentado, “Superconducting parametric amplifiers: The state of the art in Josephson parametric amplifiers,” _IEEE Microwave Magazine_ , vol. 21, no. 8, pp. 45–59, Jul 2020.
  * [8] H. Zimmer, “Parametric amplification of microwaves in superconducting Josephson tunnel junctions,” _Applied Physics Letters_ , vol. 10, no. 7, pp. 193–195, 04 1967.
  * [9] P. Russer, “Parametric amplification with Josephson junctions,” _AEÜ Archiv der Elektrischen Übertragung_ , vol. 23, no. 8, pp. 417–420, 1969.
  * [10] M. A. Castellanos-Beltran and K. W. Lehnert, “Widely tunable parametric amplifier based on a superconducting quantum interference device array resonator,” _Applied Physics Letters_ , vol. 91, no. 8, p. 083509, Aug 2007\.
  * [11] M. Esposito, A. Ranadive, L. Planat, and N. Roch, “Perspective on traveling wave microwave parametric amplifiers,” _Applied Physics Letters_ , vol. 119, no. 12, p. 120501, Sep 2021.
  * [12] K. O’Brien, C. Macklin, I. Siddiqi, and X. Zhang, “Resonant phase matching of Josephson junction traveling wave parametric amplifiers,” _Phys. Rev. Lett._ , vol. 113, p. 157001, Oct 2014.
  * [13] C. Macklin, K. O’Brien, D. Hover, M. E. Schwartz, V. Bolkhovsky, X. Zhang, W. D. Oliver, and I. Siddiqi, “A near–quantum-limited Josephson traveling-wave parametric amplifier,” _Science_ , vol. 350, no. 6258, pp. 307–310, Sep 2015.
  * [14] T. C. White, J. Y. Mutus, I.-C. Hoi, R. Barends, B. Campbell, Y. Chen, Z. Chen, B. Chiaro, A. Dunsworth, E. Jeffrey _et al._ , “Traveling wave parametric amplifier with Josephson junctions using minimal resonator phase matching,” _Applied Physics Letters_ , vol. 106, no. 24, p. 242601, Jun 2015\.
  * [15] A. L. Grimsmo and A. Blais, “Squeezing and quantum state engineering with Josephson travelling wave amplifiers,” _npj Quantum Information_ , vol. 3, p. 20, Jun 2017.
  * [16] J. Y. Qiu, A. Grimsmo, K. Peng, B. Kannan, B. Lienhard, Y. Sung, P. Krantz, V. Bolkhovsky, G. Calusine, D. Kim _et al._ , “Broadband squeezed microwaves and amplification with a Josephson travelling-wave parametric amplifier,” _Nature Physics_ , vol. 19, pp. 706–713, May 2023.
  * [17] L. Planat, A. Ranadive, R. Dassonneville, J. Puertas Martínez, S. Léger, C. Naud, O. Buisson, W. Hasch-Guichard, D. M. Basko, and N. Roch, “Photonic-crystal Josephson traveling-wave parametric amplifier,” _Phys. Rev. X_ , vol. 10, p. 021021, Apr 2020.
  * [18] A. Ranadive, M. Esposito, L. Planat, E. Bonet, C. Naud, O. Buisson, W. Guichard, and N. Roch, “Kerr reversal in Josephson meta-material and traveling wave parametric amplification,” _Scientific Reports_ , vol. 13, Apr 2022.
  * [19] A. Osman, J. Simon, A. Bengtsson, S. Kosen, P. Krantz, D. P. Lozano, M. Scigliuzzo, P. Delsing, J. Bylander, and A. Fadavi Roudsari, “Simplified Josephson-junction fabrication process for reproducibly high-performance superconducting qubits,” _Applied Physics Letters_ , vol. 118, no. 6, p. 064002, Feb 2021.
  * [20] A. Pishchimova, N. Smirnov, D. Ezenkova, E. Krivko, E. Zikiy, D. Moskalev, A. Ivanov, N. Korshakov, and I. Rodionov, “Improving Josephson junction reproducibility for superconducting quantum circuits: Junction area fluctuation,” _Scientific Reports_ , vol. 13, Apr 2023.
  * [21] J. Bylander, “Reproducible Josephson junctions are important for scalable quantum technologies,” _Superconductor Science and Technology_ , vol. 33, no. 11, p. 110501, Sep 2020.
  * [22] D. C. Feng, M. Vahidpour, Y. Mohan, N. Sharac, T. Whyland, S. Stanwyck, G. Ramachandran, and M. Selvanayagam, “Design and measurement of a Josephson traveling wave parametric amplifier fabricated in a superconducting qubit process,” in _2020 IEEE/MTT-S International Microwave Symposium (IMS)_ , 2020, pp. 940–943.
  * [23] O. Yaakobi, L. Friedland, C. Macklin, and I. Siddiqi, “Parametric amplification in Josephson junction embedded transmission lines,” _Phys. Rev. B_ , vol. 87, p. 144301, Apr 2013.
  * [24] A. B. Zorin, “Josephson traveling-wave parametric amplifier with three-wave mixing,” _Phys. Rev. Appl._ , vol. 6, p. 034006, Sep 2016.
  * [25] Y. Yuan, M. Haider, J. A. Russer, P. Russer, and C. Jirauschek, “Circuit quantum electrodynamic model of dissipative-dispersive Josephson traveling-wave parametric amplifiers,” _Phys. Rev. A_ , vol. 107, p. 022612, Feb 2023.
  * [26] M. T. Bell and A. Samolov, “Traveling-wave parametric amplifier based on a chain of coupled asymmetric squids,” _Phys. Rev. Appl._ , vol. 4, p. 024014, Aug 2015.
  * [27] T. Dixon, J. Dunstan, G. Long, J. Williams, P. Meeson, and C. Shelly, “Capturing complex behavior in Josephson traveling-wave parametric amplifiers,” _Phys. Rev. Appl._ , vol. 14, p. 034058, Sep 2020.
  * [28] S. O. Peatáin, T. Dixon, P. J. Meeson, J. M. Williams, S. Kafanov, and Y. A. Pashkin, “Simulating the effects of fabrication tolerance on the performance of Josephson travelling wave parametric amplifiers,” _Superconductor Science and Technology_ , vol. 36, no. 4, p. 045017, Mar 2023.
  * [29] C. Guarcello, G. Avallone, C. Barone, M. Borghesi, S. Capelli, G. Carapella, A. P. Caricato, I. Carusotto, A. Cian, D. Di Gioacchino _et al._ , “Modeling of Josephson traveling wave parametric amplifiers,” _IEEE Transactions on Applied Superconductivity_ , vol. 33, no. 1, pp. 1–7, Jan 2023\.
  * [30] S. T. Elkin, M. Haider, and T. E. Roth, “Hybrid numerical modeling method for Josephson traveling wave parametric amplifiers,” in _2024 IEEE International Symposium on Antennas and Propagation and USNC-URSI Radio Science Meeting (AP-S/URSI)_. (accepted for publication), 2024.
  * [31] M. Stone and P. Goldbart, _Mathematics for Physics: A Guided Tour for Graduate Students_. Cambridge University Press, 2009.
  * [32] T. E. Roth and S. T. Elkin, “Maxwell-Schrödinger modeling of a superconducting qubit coupled to a transmission line network,” _IEEE Journal on Multiscale and Multiphysics Computational Techniques_ , vol. 9, pp. 61–74, Jan 2024.
  * [33] S. M. Girvin, “Circuit QED: Superconducting qubits coupled to microwave photons,” in _Quantum Machines: Measurement and Control of Engineered Quantum Systems_. Oxford University Press, 2014, pp. 113–256.
  * [34] J.-M. Jin, _Theory and Computation of Electromagnetic Fields_. John Wiley & Sons, Ltd, 2010.
  * [35] R. A. Chilton and R. Lee, “The discrete origin of FETD-Newmark late time instability, and a correction scheme,” _Journal of Computational Physics_ , vol. 224, no. 2, pp. 1293–1306, Jun 2007.
  * [36] M. Tinkham, _Introduction to Superconductivity: Second Edition_ , ser. Dover Books on Physics. Dover Publications, 2004.
  * [37] A. Zorin, “Flux-driven Josephson traveling-wave parametric amplifier,” _Phys. Rev. Appl._ , vol. 12, p. 044051, Oct 2019.
  * [38] M. Haider, “A black-box quantum model for superconducting traveling-wave parametric amplifiers,” _IEEE Transactions on Microwave Theory and Techniques_ , pp. 1–15, Jan 2024.

|  Samuel T. Elkin (S’22) received the B.S degree in electrical engineering
from Purdue University, West Lafayette, IN, USA, in 2022. From 2022 to 2023,
he worked at Indesign, LLC., Indianapolis, IN, USA, as an electrical engineer.
He is currently pursuing a Ph.D. degree in electrical and computer engineering
at Purdue University. His current research interests are in multiphysics
modeling of superconducting circuit quantum devices.
---|---
|  Michael Haider (Member, IEEE) received the B.Sc., M.Sc., and Ph.D. (summa
cum laude) degrees in electrical and computer engineering from the Technical
University of Munich (TUM), Garching, Germany, in 2014, 2016, and 2019,
respectively. From 2016 to 2019, he was a Research Associate with the
Institute of Nanoelectronics, TUM, where his research focused on quantum
circuit theory, wireless power transfer, and stochastic electromagnetic
fields. In 2019, he joined the Computational Photonics Group, TUM, as a Post-
Doctoral Researcher, where he became a permanent staff Scientist in October
2023. He has coauthored over 100 scientific papers in refereed journals and
conference proceedings. His current research interests include the modeling of
charge carrier transport in quantum cascade lasers, the dynamical modeling of
quantum devices within a generalized Maxwell-Bloch framework, and quantum
circuit theoretical models applied to Josephson traveling-wave parametric
amplifiers. Dr. Haider received the 2020 VDE ITG Dissertationspreis for his
outstanding Ph.D. thesis on Investigations of Stochastic Electromagnetic
Fields.
---|---
|  Thomas E. Roth (Member, IEEE) received the B.S. degrees in electrical
engineering and computer engineering from the Missouri University of Science
and Technology, Rolla, MO, USA, in 2015, and the M.S. and Ph.D. degrees in
electrical and computer engineering from the University of Illinois at Urbana-
Champaign, Champaign, IL, USA, in 2017 and 2020, respectively. He is currently
an Assistant Professor with the Elmore Family School of Electrical and
Computer Engineering, Purdue University, West Lafayette, IN, USA. Prior to
joining Purdue, he was a Senior Member of the Technical Staff with Sandia
National Laboratories, Albuquerque, NM, USA. He has received five best paper
or young scientist awards from conferences, and the Ruth and Joel Spira
Outstanding Teacher Award from Purdue University. His current research
interests include multiscale and multiphysics computational electromagnetics
techniques, particularly for rigorously analyzing quantum information
processing devices.
---|---


# CAMIRADA: Cancer microRNA association discovery algorithm, a case study on
breast cancer

Sepideh Shamsizadeh<EMAIL_ADDRESS>Sama Goliaei<EMAIL_ADDRESS>Zahra Razaghi Moghadam<EMAIL_ADDRESS>Faculty of New Sciences and
Technologies, University of Tehran, Tehran, Iran  Max Planck Institute of
Molecular Plant Physiology, Posdam, German

###### Abstract

In recent studies, non-coding protein RNAs have been identified as microRNA
that can be used as biomarkers for early diagnosis and treatment of cancer,
that decrease mortality in cancer. A microRNA may target hundreds or thousands
of genes and a gene may regulate several microRNAs, so determining which
microRNA is associated with which cancer is a big challenge. Many
computational methods have been performed to detect micoRNAs association with
cancer, but more effort is needed with higher accuracy. Increasing research
has shown that relationship between microRNAs and TFs play a significant role
in the diagnosis of cancer. Therefore, we developed a new computational
framework (CAMIRADA) to identify cancer-related microRNAs based on the
relationship between microRNAs and disease genes (DG) in the protein network,
the functional relationships between microRNAs and Transcription Factors (TF)
on the co-expression network, and the relationship between microRNAs and the
Differential Expression Gene (DEG) on co-expression network. The CAMIRADA was
applied to assess breast cancer data from two HMDD and miR2Disease databases.
In this study, the AUC for the 65 microRNAs of the top of the list was 0.95,
which was more accurate than the similar methods used to detect microRNAs
associated with the cancer artery.

###### keywords:

microRNAs , Transcription factors , Differentially expressed gene , Diseases
gene , Co-expression network

††journal: Journal of Biomedical Informatics

##  1 Introduction

According to a study,breast cancer is one of the most common causes of women’s
deaths in cancer [1]. The molecular traits of the primary tumors play an
important role in timely diagnosis and treatment of the next stages of breast
cancer. Therefore, we need methods that can detect cancer at an early stage in
order to help clinicians and patients to be treated with biological markers
[1]. In recent years, many studies have been conducted to clarify the
mechanisms that make cancer progress and develop.

Although many genes that cause cancer and stop it are known by researchers, it
is still necessary to identify cancer pathways. Therefore, one of the most
important biological targets of cancer is the detection of genes related to
cancers [2, 3].

MicroRNAs are short non-coding RNAs with an approximate length of 22
nucleotides that are involved in post-transcriptional regulation and that are
major regulators of gene expression [4].In 2002, the first link between cancer
and microRNAs was identified [5]. MicroRNAs target about 60% of human genes
that be involved in a wide range of biological processes and disease including
cell division, proliferation, differentiation, and apoptosis [6]. More than
50% of human microRNAs are located in cancer-associated genomic regions.
MicroRNAs could play roles as oncogenes or tumor suppressor genes.
Dysregulation of microRNAs expression has been shown to have impacts on human
diseases [6].

However, identifying the related microRNAs with existing experimental tests
may be difficult and time consuming. In addition, many researchers are also
faced with limited knowledge about microRNAs. Therefore, a number of
computational approaches have been recently developed to identify microRNAs
associated with the disease. Statistical methods, machine learning, and
network-based methods are some of the methods proposed to predict cancer-
related microRNAs.

Zhao et al. developed a framework for obtaining the relationship between
cancer and microRNAs through their target gene expression profiles without
requiring neither microRNA expression data or the matched gene and microRNAs
expression data. For each microRNA, target genes that are likely to be related
to cancer are determined, and subsequently identified adverse pathways
associated with cancer. Then, microRNAs are ranked according to these
pathways, because high-ranking microRNAs are more likely to be linked to
cancer [7].

Xue et al. used the Gaussian mixture models, they identified patterns of gene
expression for healthy and patient samples, and then, with Fisher’s exact
test, deduced the regulator relationship between TF and genes. A minimum
description length for the pruning of the network, using the (MDL) principle,
the relationships between TFs and microRNAs have been achieved [8].

Tseng et al. provided a method to identify active oncomirs and their potential
functions in gastric cancer progression. The microRNA and mRNA’s expression
profiles with the human protein interaction network (PIN) are integrated to
show microRNAs-regulated PIN in specific biological conditions. The microRNAs’
potential functions were identified by functional enrichment analysis and the
activities of microRNAs-regulated PINs were evaluated by the co-expression of
protein-protein interactions (PPIs) [9].

Le constructed microRNAs functional similarity networks based on shared
targets of microRNAs, and then it integrated them with a microRNA functional
synergistic network. After analyzing topological properties of these networks,
it introduced five network-based ranking methods (RWR, PRINCE (Prioritization
and Complex Elucidation), which was proposed for disease gene prediction;
PageRank with Priors (PRP) and K-Step Markov (KSM), which were used for
studying web networks; and a neighborhood-based algorithm.) in prioritizing
candidate microRNAs to predict novel disease-related microRNAs based on the
constructed microRNAs functional similarity networks [6].

In this study, we presented CAMIRADA framework for the identification of
microRNAs cancer associations by taking advantage of the useful links between
microRNAs targets and disease genes in protein-protein networks (PPIs), the
functional connections between microRNAs targets and Transcription Factors(TF)
in co-expression network and the functional relationships between microRNAs
targets and Differentially Expressed Gene(DEG) in co-expression network.

##  2 Methods

###  2.1 Data Sets

We firstly collected the microRNAs target genes predicted with different
tools, including PicTar [10], miRanda (version 3.0) [11], TargetScan (release
6.2) [12], miRBase (version 16) [13] and mirTarbase [14]. Specifically, for a
microRNA, we consider target genes to be obtained at least by three tools.
Note that some microRNAs are not considered if, after filtering, its target
genes are less than three. In total, we obtained 42832 targeting pairs that
involved 825 microRNAs and 8334 target genes.

The disease-gene association data were obtained from DisGeNET [15] and The
Candidate Cancer Gene Database (CCGD) [16], which involved 375 genes related
to breast cancer. The Transcription factors(TF) were retrieved from ChIPBase
[17], that included 13890 genes. The names of all these genes mapped to Entrez
gene id.

###  2.2 Human Protein-Protein Interaction (PPI) Data And Co-Expression
Network

PPI data for humans are derived from the Human Protein Reference Database
(HPRD Release 9), which contains human-protein annotations based on empirical
evidence from published reports [18]. We changed the name of the gene to the
Entrez gene id, and then we obtained the maximum components of the entire
network, which contains 9028 genes and 35865 interactions. It is noteworthy
that PPI data in HPRD are the most common protein isoforms, mainly due to lack
of experimental data [12].

We downloaded GSE31192 [19] data from the GEO [20] database to build a co-
expression network and obtain DEGs. The RNA extracted from breast tissue with
microarrays Affymetrix Human Genome U133 Plus 2.0 has been investigated. These
gene expression data are pre-processed with the Robust Multichip Average
algorithm(RMA).

Then we find out the Differentially Expressed Genes (DEGs) between normal and
abnormal samples with samr package in R. We selected genes, which their
$p-value$ were under 0.05 and their $nfold$ was 0.1 . The co-expression
network was constructed using the WGCNA [21] package in R software. To choose
an appropriate cutoff to include a percentage of the highest correlations, the
approximate scale-free topology criterion was applied. For determining the
optimal parameter, the function pickSoftThresholdthe in WGCNA package was
used. After selecting the best threshold (0.85), adjacency network for co-
expression network was calculated.

The co-expression network included 21179 genes, so for sampling, we considered
DEGs and all targets of microRNAs as seed, we mapped them onto the co-
expression network, and after that, we applied RWR (Random Walk With Restart)
algorithm for ranking the genes of network [18]. Then we selected top 10000
genes and made the network base on their relationships.

###  2.3 CAMIRADA algorithm

####  2.3.1 Disease Genes-microRNAs

We mapped targets of a microRNA and disease genes onto PPI network. After
that, we considered disease genes as seeds and applied RWR [18] to rank genes
of PPI network. Then by using gene set enrichment analysis (GSEA) we derived
the ranked gene list, we used ES1 (enrichment score) using the following
formula [18].

$\text{ES1}=\max_{1\leq i\leq n}\left(\sum_{\begin{subarray}{c}g_{j}\in TG,\\\
j\leq
i\end{subarray}}\sqrt{\frac{N-n_{1}}{n_{1}}}-\sum_{\begin{subarray}{c}g_{j}\notin
TG,\\\ j\leq i\end{subarray}}\sqrt{\frac{N-n_{1}}{n_{1}}}\right)$ (1)

Where $N$ is the number of all genes in PPI network and $n_{1}$ is the number
of target genes of one microRNA that we showed with TG. For each microRNA we
calculate ES1. We computed a statistical sum from the beginning of the ranked
list of PPI’s genes. In this way, by moving down the list if the gene was in
TG, we added to statistical sum, and if that gene did not exist in the TG, we
reduced from statistical sum. The RWR algorithm with each microRNA’s target
genes as seeds was applied to compute ES2 for the same pairing of microRNAs
disease referred to above [18].

$\text{ES2}=\max_{1\leq i\leq n}\left(\sum_{\begin{subarray}{c}g_{j}\in DG,\\\
j\leq
i\end{subarray}}\sqrt{\frac{N-n_{1}}{n_{1}}}-\sum_{\begin{subarray}{c}g_{j}\notin
DG,\\\ j\leq i\end{subarray}}\sqrt{\frac{N-n_{1}}{n_{1}}}\right)$ (2)

Which $n_{2}$ is the number of disease genes. For all pairs of microRNAs and
disease genes we did the above procedures. For computing ES1 we considered
disease genes as seed for RWR algorithm and we considered targets of one
microRNA as seed for RWR to compute ES2. Then we calculated ES for each pair
of microRNAs-disease by using the following formula [18]

$\text{ES}=\beta\ \text{ES1}+(1-\beta)\ \text{ES2}$ (3)

####  2.3.2 DEGs-microRNAs

We mapped targets of each microRNA and DEG onto co-expression network. For all
pairs of microRNAs and DEGs we calculated ES. For computing ES1 we considered
DEGs as seed for RWR algorithm and we considered targets of one microRNA as
seed for RWR to compute ES2. Then we calculated ES for each pair of microRNAs-
DEG Where $N$ is the number of all genes in co-expression network and $n_{1}$
is the number of target genes of one microRNA. For each microRNA we calculate
ES1 1. Which $n_{2}$ in the ES2 2 is the number of DEGs.

Figure 1: An overview of our framework to predict new diseases associated with
microRNAs.

####  2.3.3 TFs-microRNAs

We did all above producers on TF and targets of microRNAs. First in order, we
mapped targets of all microRNAs and TFs onto co-expression network and for
each microRNA we separately saved it’s targets in a list. For all pairs of
microRNAs and TFs we calculated ES 3. For computing ES1 we considered TFs as
seed for RWR algorithm and we considered targets of one microRNA as seed for
RWR to compute ES2. Then we calculated ES for each pair of microRNAs-TFs Where
$N$ is the number of all genes in co-expression network and $n_{1}$ is the
number of target genes of one microRNA. For each microRNA we calculate ES1 1.
Which $n_{2}$ in the ES2 2 is the number of TFs. A review of our frame work
are showed in Figure 1.

###  2.4 Random Network And P-value

To measure the importance of the relationship between microRNAs and disease,
we used the p-value. To calculate the p-value, we first created 1000 random
networks of the PPI network and 1000 random networks of the co-expression
network. Random networks were made in such a way that the number of input and
output edges from each node with the number of input and output edges from the
same node in the primary network are equal [18]. Of course, it should be noted
that both the PPI network and the co-expression network are non-directional
networks in this study. The p-value is calculated based on [18].

##  3 Results And Performance

To evaluate the performance of our algorithm to identify microRNAs-disease
associations, we plotted the receiver operating characteristic (ROC) curves
and computed the area under the curve (AUC).

To compute the AUC, we first need to calculate the sensitivity and specificity
based on the following formulas, and the sensitivity is defined as the correct
detection rate of the positive group on the $y$ axis and the specificity, or
the wrong detection rate of the negative category on the $x$ axis. A ROC curve
shows a relative compromise between profits and costs [6].

$1-\text{specificity}=\frac{FP}{TN+FP}$ (4)
$1-\text{Sensitivity}=\frac{TP}{FN+TP}$ (5)

For computing of these we needed to define the set of $TP$، $TN$، $FP$ and
$FN$. At first we should define positive set and negative set. For positive
set we considered the known cancer related microRNAs were obtained from
miR2Disease [22] and HMDD [23] databases.

Nowadays, it is very difficult or impossible to collect data for non-cancerous
microRNA [24]. For a negative data set, we chose microRNAs that exhibited the
lowest fold change values as negative controls by analyzing the corresponding
expression profile of the breast cancer. We also used the same number of
negative controls as positive ones. We downloaded microRNAs expression profile
GSE45666 [24] from GEO.

We computed three lists of p-value for microRNAs according to our method. We
arranged lists of p-value by descending order and we calculated AUC for each
list to identify the importance of each factor (TF, disease and DEG) for
determining which microRNA is associated with breast cancer.

As the Figure 2 shows, the lowest AUC is for disease-microRNAs and the maximum
AUC is for CAMIRADA approach. So, this shows that DEGs and TFs influence on
identification of which microRNAs are related to breast cancers. For using the
influence of TFs and DEGs on identification of microRNAs, we combined p-vlues
with some tests in R.

Figure 2: Computing AUC for each three list of p-value for DEG-microRNAs,
Disease-microRNAs and TF-microRNAs. And we calculated AUC for CAMIRADA.

The tests included Wilkinsonp, sumlog, sumz, logitp, meanp and sump from metap
[25] package in R. Then we scored microRNAs according to our model. For each
p-value as threshold, we computed $TP$, $TN$, $FP$ and $FN$ by using the
following definition:

  * 1.

TP: The number of microRNAs are in positive set and their scores are below the
threshold.

  * 2.

FP: The number of microRNAs are in negative set or they are not in positive
set and their scores are below the threshold.

  * 3.

TN: The number of microRNAs are in negative set or they are not in positive
and their scores are above the threshold.

  * 4.

FN: The number of microRNAs are in positive set and their scores are above the
threshold.

Figure 3: The AUC for all of the p-value combining tests.

To select the best test for combining the p-value, we combined the three lists
of p-values by all tests and then we calculated AUC for each tests results are
shown in Figure 3. As shown in Figure 3, the best result is for sumlog test by
82%, so we choose this test for combining the three lists of p-values.

For more comparing the tests to combine the p-value we combined the three
lists of p-values and we calculated AUC for 65, 130, and 260 number of the
microRNAs from the top of the ranked list of microRNAs and the result of the
comparison of tests are shown by bar chart in Figure 4.

Figure 4: Bar chart for comparing the tests to combine the p-value.

According the results that achieved from tests of p-value combining, sumlog
choose to combine p-value. So, we combined three ranked list of p-value with
sumlog and then arranged the list of combined p-vlaues in ascending order. The
lower the p-value for a microRNA, the greater the probability that this
microRNA is associated with the cancer.

In Table 1 the 10 top of the list of the ranked microRNAs with CAMIRADA and
the role of microRNAs in breast cancer are shown.

Table 1: The 10 top of the list of the ranked microRNAs. microRNAs | p-value | role of microRNAs in breast cancer
---|---|---
hsa-miR-93 [26] | 0.000679 | Tumor suppressor
hsa-miR-526b [27] | 0.000702 | Oncogene
hsa-miR-17 [28] | 0.000951 | Oncogene
hsa-miR-20b [29] | 0.00125 | Oncogene
hsa-miR-199a [30] | 0.001486 | Tumor suppressor
hsa-miR-20a [31] | 0.001529 | Oncogene
hsa-miR-149 [32] | 0.0020011 | Oncogene
hsa-miR-199b [33] | 0.00265 | Biomarker
hsa-miR-519d [38] | 0.002863 | Oncogene
hsa-miR-502 [34] | 0.036284 | Tumor suppressor

##  4 Conclusion

The identification of novel cancer-associated microRNAs is important for
diagnosing cancer in primary steps for treatment the cancer. In this study, we
presented a new approach (CAMIRADA) to identify microRNAs related to breast
cancer. CAMIRADA find the relationship between microRNAs and cancer and it
ranked microRNAs base on their relation with cancer, by using the microRNAs-
disease genes associated, the DEGs-microRNAs related, and the TFs-microRNAs
relation. To evaluated CAMIRADA algorithm we calculated AUC which was 0.95 for
65 microRNAs in top of the ranked list of microRNAs. CAMIRADA showed that TFs
and DEGs have important roles for identification of which microRNA is related
to cancer, so for showing the influence of them we combined three ranked lists
of p-values by sumlog test.

By using CAMIRADA, we predicted eight novel microRNAs including hsa-miR-93,
hsa-miR-526b, hsa-miR-20b, hsa-miR-199a, hsa-miR-6884, hsa-miR-199b, and hsa-
miR-519d associated to breast cancer which are not yet recorded in the
disease-microRNA association HMDD and miR2Disease. For future studies, our
methods can be used to identify which microRNAs are related which cancers, we
can make a network for each microRNA and associated cancers, then analysis the
network and extract information to identify related between cancers.

## References

## References

  * Yang et al. [2015] Yang Yang, Yiqiao Xing, Chaoqun Liang, Liya Hu, Fei Xu, and Yuan Chen.  Crucial micrornas and genes of human primary breast cancer explored by microrna-mrna integrated analysis.  _Tumor Biology_ , 36(7):5571–5579, 2015.
  * Farahmand et al. [2016] S Farahmand, S Goliaei, N Ansari-Pour, and Z Razaghi-Moghadam.  Gta: a game theoretic approach to identifying cancer subnetwork markers.  _Molecular BioSystems_ , 12(3):818–825, 2016\.
  * Razaghi-Moghadam et al. [2016] Zahra Razaghi-Moghadam, Razieh Abdollahi, Sama Goliaei, and Morteza Ebrahimi.  Hybridranker: Integrating network topology and biomedical knowledge to prioritize cancer candidate genes.  _Journal of biomedical informatics_ , 64:139–146, 2016\.
  * Guruceaga and Segura [2014] Elizabeth Guruceaga and Victor Segura.  Functional interpretation of microrna–mrna association in biological systems using r.  _Computers in biology and medicine_ , 44:124–131, 2014\.
  * Adalı and Şekeroğlu [2012] Terin Adalı and Boran Şekeroğlu.  Analysis of micrornas by neural network for early detection of cancer.  _Procedia Technology_ , 1:449–452, 2012.
  * Le [2015] Duc-Hau Le.  Network-based ranking methods for prediction of novel disease associated micrornas.  _Computational Biology and Chemistry_ , 58:139–148, 2015\.
  * Zhao et al. [2014] Xing-Ming Zhao, Ke-Qin Liu, Guanghui Zhu, Feng He, Béatrice Duval, Jean-Michel Richer, De-Shuang Huang, Chang-Jun Jiang, Jin-Kao Hao, and Luonan Chen.  Identifying cancer-related micrornas based on gene expression data.  _Bioinformatics_ , 31(8):1226–1234, 2014.
  * Xue et al. [2017] Mengzhu Xue, Haiyue Liu, Liwen Zhang, Hongyuan Chang, Yuwei Liu, Shaowei Du, Yingqun Yang, and Peng Wang.  Computational identification of mutually exclusive transcriptional drivers dysregulating metastatic micrornas in prostate cancer.  _Nature communications_ , 8:14917, 2017.
  * Tseng et al. [2011] Chien-Wei Tseng, Chen-Ching Lin, Chiung-Nien Chen, Hsuan-Cheng Huang, and Hsueh-Fen Juan.  Integrative network analysis reveals active micrornas and their functions in gastric cancer.  _BMC systems biology_ , 5(1):99, 2011.
  * Krek et al. [2005] Azra Krek, Dominic Grün, Matthew N Poy, Rachel Wolf, Lauren Rosenberg, Eric J Epstein, Philip MacMenamin, Isabelle Da Piedade, Kristin C Gunsalus, Markus Stoffel, et al.  Combinatorial microrna target predictions.  _Nature genetics_ , 37(5):495, 2005.
  * John et al. [2004] Bino John, Anton J Enright, Alexei Aravin, Thomas Tuschl, Chris Sander, and Debora S Marks.  Human microrna targets.  _PLoS biology_ , 2(11):e363, 2004.
  * Lewis et al. [2005] Benjamin P Lewis, Christopher B Burge, and David P Bartel.  Conserved seed pairing, often flanked by adenosines, indicates that thousands of human genes are microrna targets.  _cell_ , 120(1):15–20, 2005.
  * Griffiths-Jones et al. [2006] Sam Griffiths-Jones, Russell J Grocock, Stijn Van Dongen, Alex Bateman, and Anton J Enright.  mirbase: microrna sequences, targets and gene nomenclature.  _Nucleic acids research_ , 34(suppl_1):D140–D144, 2006.
  * Piñero et al. [2016] Janet Piñero, Àlex Bravo, Núria Queralt-Rosinach, Alba Gutiérrez-Sacristán, Jordi Deu-Pons, Emilio Centeno, Javier García-García, Ferran Sanz, and Laura I Furlong.  Disgenet: a comprehensive platform integrating information on human disease-associated genes and variants.  _Nucleic acids research_ , page gkw943, 2016.
  * Abbott et al. [2014] Kenneth L Abbott, Erik T Nyre, Juan Abrahante, Yen-Yi Ho, Rachel Isaksson Vogel, and Timothy K Starr.  The candidate cancer gene database: a database of cancer driver genes from forward genetic screens in mice.  _Nucleic acids research_ , 43(D1):D844–D848, 2014\.
  * Yang et al. [2012] Jian-Hua Yang, Jun-Hao Li, Shan Jiang, Hui Zhou, and Liang-Hu Qu.  Chipbase: a database for decoding the transcriptional regulation of long non-coding rna and microrna genes from chip-seq data.  _Nucleic acids research_ , 41(D1):D177–D187, 2012\.
  * Chou et al. [2015] Chih-Hung Chou, Nai-Wen Chang, Sirjana Shrestha, Sheng-Da Hsu, Yu-Ling Lin, Wei-Hsiang Lee, Chi-Dung Yang, Hsiao-Chin Hong, Ting-Yen Wei, Siang-Jyun Tu, et al.  mirtarbase 2016: updates to the experimentally validated mirna-target interactions database.  _Nucleic acids research_ , 44(D1):D239–D247, 2015\.
  * Shi et al. [2013] Hongbo Shi, Juan Xu, Guangde Zhang, Liangde Xu, Chunquan Li, Li Wang, Zheng Zhao, Wei Jiang, Zheng Guo, and Xia Li.  Walking the interactome to identify human mirna-disease associations through the functional link between mirna targets and disease genes.  _BMC systems biology_ , 7(1):101, 2013.
  * Harvell et al. [2013] Djuana ME Harvell, Jihye Kim, Jenean O’Brien, Aik-Choon Tan, Virginia F Borges, Pepper Schedin, Britta M Jacobsen, and Kathryn B Horwitz.  Genomic signatures of pregnancy-associated breast cancer epithelia and stroma and their regulation by estrogens and progesterone.  _Hormones and Cancer_ , 4(3):140–153, 2013.
  * Barrett et al. [2012] Tanya Barrett, Stephen E Wilhite, Pierre Ledoux, Carlos Evangelista, Irene F Kim, Maxim Tomashevsky, Kimberly A Marshall, Katherine H Phillippy, Patti M Sherman, Michelle Holko, et al.  Ncbi geo: archive for functional genomics data sets—update.  _Nucleic acids research_ , 41(D1):D991–D995, 2012\.
  * Langfelder and Horvath [2008] Peter Langfelder and Steve Horvath.  Wgcna: an r package for weighted correlation network analysis.  _BMC bioinformatics_ , 9(1):559, 2008.
  * Jiang et al. [2008] Qinghua Jiang, Yadong Wang, Yangyang Hao, Liran Juan, Mingxiang Teng, Xinjun Zhang, Meimei Li, Guohua Wang, and Yunlong Liu.  mir2disease: a manually curated database for microrna deregulation in human disease.  _Nucleic acids research_ , 37(suppl_1):D98–D104, 2008.
  * Li et al. [2013] Yang Li, Chengxiang Qiu, Jian Tu, Bin Geng, Jichun Yang, Tianzi Jiang, and Qinghua Cui.  Hmdd v2. 0: a database for experimentally supported human microrna and disease associations.  _Nucleic acids research_ , 42(D1):D1070–D1074, 2013.
  * Lee et al. [2013] Chia-Hsien Lee, Wen-Hong Kuo, Chen-Ching Lin, Yen-Jen Oyang, Hsuan-Cheng Huang, and Hsueh-Fen Juan.  Microrna-regulated protein-protein interaction networks and their functions in breast cancer.  _International journal of molecular sciences_ , 14(6):11560–11606, 2013.
  * [25] Metap, a program to combine p values.  URL http://igm.cumc.columbia.edu/MetaP/metap.php.
  * Shyamasundar et al. [2016] Sukanya Shyamasundar, Jia Pei Lim, and Boon Huat Bay.  mir-93 inhibits the invasive potential of triple-negative breast cancer cells in vitro via protein kinase wnk1.  _International journal of oncology_ , 49(6):2629–2636, 2016.
  * Landman [2014] Erin O Landman.  The role of mir-526b in cox-2 mediated human breast cancer progression and induction of stem-like phenotype via ep4 receptor signaling.  2014\.
  * Hossain et al. [2006] Anwar Hossain, Macus T Kuo, and Grady F Saunders.  Mir-17-5p regulates breast cancer cell proliferation by inhibiting translation of aib1 mrna.  _Molecular and cellular biology_ , 26(21):8191–8201, 2006.
  * Li et al. [2012] Jian-Yi Li, Yang Zhang, Wen-Hai Zhang, Shi Jia, Ye Kang, and Xiao-Yu Zhu.  Differential distribution of mir-20a and mir-20b may underly metastatic heterogeneity of breast cancers.  _Asian Pacific Journal of Cancer Prevention_ , 13(5):1901–1906, 2012.
  * Chen et al. [2016] Jiawei Chen, Vivian Y Shin, Man T Siu, John CW Ho, Isabella Cheuk, and Ava Kwong.  mir-199a-5p confers tumor-suppressive role in triple-negative breast cancer.  _BMC cancer_ , 16(1):887, 2016.
  * Li et al. [2015] Pu Li, Jing-Xuan Shan, Xue-Hua Chen, Di Zhang, Li-Ping Su, Xiu-Ying Huang, Bei-Qin Yu, Qiao-Ming Zhi, Cheng-Long Li, Ya-Qing Wang, et al.  Epigenetic silencing of microrna-149 in cancer-associated fibroblasts mediates prostaglandin e2/interleukin-6 signaling in the tumor microenvironment.  _Cell research_ , 25(5):588, 2015.
  * Fang et al. [2016] Cheng Fang, Fu-Bing Wang, Yirong Li, and Xian-Tao Zeng.  Down-regulation of mir-199b-5p is correlated with poor prognosis for breast cancer patients.  _Biomedicine & Pharmacotherapy_, 84:1189–1193, 2016.
  * Lehmann et al. [2010] Ulrich Lehmann, Thomas Streichert, Benjamin Otto, Cord Albat, Britta Hasemeier, Henriette Christgen, Elisa Schipper, Ursula Hille, Hans H Kreipe, and Florian Länger.  Identification of differentially expressed micrornas in human male breast cancer.  _BMC cancer_ , 10(1):109, 2010.
  * Cheng et al. [2012] Chun-Wen Cheng, Hsiao-Wei Wang, Chia-Wei Chang, Hou-Wei Chu, Cheng-You Chen, Jyh-Cherng Yu, Jui-I Chao, Huei-Fang Liu, Shian-ling Ding, and Chen-Yang Shen.  Microrna-30a inhibits cell migration and invasion by downregulating vimentin expression and is a potential prognostic marker in breast cancer.  _Breast cancer research and treatment_ , 134(3):1081–1093, 2012.