pnasresearcharticle Majumder Wave breaks occur in many complex excitable
systems and often lead to self-sustained re-entry. In the heart, such re-entry
causes fatal rhythm disturbances, so-called arrhythmias, the prevention of
which requires a detailed understanding of the initiation and dynamics of
these breaks. Optogenetics provides a promising tool to investigate and
control wave dynamics in cardiac tissue. Although it is very effective at
stimulating light intensities, the tools’ performance remains largely unknown
at non-stimulating intensity levels. Here we demonstrate a new mechanism of
wave break initiation in cardiac tissue, which is triggered by periodic sub-
threshold illumination. The mechanism involves "conditioning" the wavelength
of a propagating wave, which leads to the creation of new excitation windows
for re-entry. Please provide details of author contributions here. The authors
declare no competing interests. 1Rupamanjari Majumder. E-mail:
rupamanjari.majumderds.mpg.de

# A mechanism for electric turbulence in cardiac tissue with optogenetic
modification.

Rupamanjari Majumder  Max Planck Institute for Dynamics and Self-Organization,
Göttingen, Germany  Sayedeh Hussaini  Max Planck Institute for Dynamics and
Self-Organization, Göttingen, Germany  Institute for Dynamics of Complex
Systems, University of Göttingen, D-37073 Göttingen, Germany  Vladimir S.
Zykov  Max Planck Institute for Dynamics and Self-Organization, Göttingen,
Germany  Stefan Luther  Max Planck Institute for Dynamics and Self-
Organization, Göttingen, Germany  Institute for Dynamics of Complex Systems,
University of Göttingen, D-37073 Göttingen, Germany  Eberhard Bodenschatz  Max
Planck Institute for Dynamics and Self-Organization, Göttingen, Germany
Institute for Dynamics of Complex Systems, University of Göttingen, D-37073
Göttingen, Germany  Laboratory of Atomic and Solid-State Physics and Sibley
School of Mechanical and Aerospace Engineering, Cornell University, Ithaca,
New York 14853, USA

###### Abstract

Interruptions in nonlinear wave propagation, commonly referred to as wave
breaks, are typical of many complex excitable systems. In the heart they lead
to fatal rhythm disorders, the so-called arrhythmias, which are one of the
main causes of sudden death in the industrialized world. Progress in the
treatment and therapy of cardiac arrhythmias requires a detailed understanding
of the triggers and dynamics of these wave breaks. In particular, two very
important questions are: 1) What determines the potential of a wave break to
initiate re-entry? and 2) How do these breaks evolve such that the system is
able to maintain spatiotemporally chaotic electrical activity? Here we
approach these questions numerically using optogenetics in an in silico model
of human atrial tissue that has undergone chronic atrial fibrillation (cAF)
remodelling. In the lesser known sub-threshold illumination régime, we
discover a new mechanism of wave break initiation in cardiac tissue that
occurs for gentle slopes of the restitution characteristics. This mechanism
involves "conditioning" or reshaping the wave profile from front to back, such
that, removal of the external light source causes rapid recovery of cells at
the waveback, leading to the creation of vulnerable windows for sustained re-
entry in spatially extended systems.

###### keywords:

Wave break $|$ Spiral waves $|$ Optogenetics $|$ Sub-threshold illumination
$|$ Mechanism

###### doi:

www.pnas.org/cgi/doi/10.1073/pnas.XXXXXXXXXX

This manuscript was compiled on

Spiral waves occur as short- or long-lived transients in various natural
excitable systems (1, 2, 3, 4, 5, 6). In the heart they appear as abnormal
electrical waves which are the basis of fatal cardiac arrhythmias (tachycardia
and fibrillation) (7). Regardless of the type of medium that sustains these
waves, spirals exhibit some common dynamic characteristics, such as drift (8),
meander (9, 10), anchoring (11, 12, 13, 14), detachment (15, 16, 17, 18) and
breakup (19). Spiral wave breakup or turbulence is of particular interest to
cardiac researchers because it is associated with fibrillation, the most
common precursor of stroke and sudden death. In particular, the mechanisms
underlying the occurrence of spiral wave turbulence in heart tissue are
extensively investigated because of their clinical significance and possible
implications (20).

The tendency of a spiral to break up in an excitable medium is critically
determined by the medium’s excitability. Spatiotemporal modulation of the
excitability of the medium can change the dynamics of spiral waves. Depending
on the amplitude, frequency, degree of synchronization and the spatial scale
of modulation of excitability, phenomena such as drift (21, 22), deformation
(23), block (23), meander (24, 22), breakup (23, 25) and suppression of spiral
waves (26) are known to occur. In experiments, however, excitability is
extremely difficult to control in real time and in a reversible manner. Thanks
to optogenetics, this has recently become possible to a considerable extent.
Optogenetics has revolutionized cardiac research by successfully demonstrating
supposedly pain-free real-time spatio-temporal control of spiral wave dynamics
in cardiac tissue (27, 28, 29, 30, 31, 32). With this tool one can regulate
the flow of current across cell membranes to force these cell to trigger or
suppress action potentials (33). Although there is an extensive literature
cataloguing the advantages and disadvantages of using optogenetics at supra-
threshold light intensities (i.e. intensities capable of stimulating action
potentials in individual cells or inducing waves in extended media), very
little is known about the functionality and properties of the tool at the
‘sub-threshold’ level. While some studies show that supra-threshold
disturbances of spiral waves can cause wave breaks (34), the corresponding
response at sub-threshold intensities is not known.

An interesting study by Park et al. (35) demonstrates the possibility to
achieve prolongation of the action potential duration (APD) in optogenetically
modified cardiomyocytes by applying optical stimuli in the subthreshold
illumination regime. In particular, they show that optical pulses can prolong
APD in a graded manner, thereby allowing one to ‘scuplt’ the action potential
(AP) morphology in general. In another study by Burton et al. (36), the
authors demonstrate that application of a subthreshold light pulse can lead to
an increase in conduction velocity of a propagating wave. However, this
increase only occurs in the speed of the first wave in a wave train. The
reason may be attributed to the current-voltage characteristics of the
Channelrhodopsin current, which displays an initial spike, followed by
saturation at a decreased level, when the optogenetically-mofidied cell is
continuously exposed to light.

In this study we show a new mechanism of wave break initiation in simulated,
optogenetically modified human atrial tissue in the presence of sub-threshold
illumination. In particular, we observe that the global periodic pulsing of
the domain with light of high sub-threshold intensity can lead to the
development of wave breaks by "conditioning” of the wavelength, i.e.
modulating the distribution of excitation along one arm of the mother spiral,
thus creating an excitation window for wave break and re-entry.

##  1 Results

Under the influence of a constant, uniform global illumination, we found that
an increase in light intensity (LI) led to an increase in the conduction
velocity (CV) of the propagating wave in one dimension (1D) and to a
progressive increase in the core size of the meandering spiral in two
dimensions (2D), while maintaining the shape of the tip trajectory (Fig. 1A,B
and D). The dominant frequencies of rotation of the spiral decreased with
increase in LI as illustrated in Fig. 1C. However, replacement of the time-
independent illumination with a periodic perturbation led to the emergence of
interesting dynamical effects. In particular,

Figure 1: Influence of a constant, uniform, global sub-threshold illumination
with light intensity (LI) on spiral wave dynamics in 2D. A) Tip trajectory of
the spiral at different LIs. B) Increase of the core diameter ($d_{core}$) of
the spiral with increase of LI. Here the core is defined as the outer circle
that just encloses the tip trajectory. C) Decrease in the dominant (first and
second fundamental) frequencies of the spiral with increase in LI. D)
Trajectory of the spiral wave at constant, uniform, global illumination below
the threshold with LI = 0.75 mW/mm2. E) Break up of the spiral wave at the
same LI, when the light is applied periodically with a frequency of 2.1Hz.
Green arrows point to the lines of the conduction block (shown in bold white
lines). The faint grey lines trace the past dynamic history of the tips of the
various spirals that have formed and terminated over time. Figure 2: A) APD90
and B) CV restitution curves in the absence of light (black), in the presence
of uniform global illumination with low sub-threshold light intensity
(LI=0.25mW/mm2, red) and in the presence of uniform global illumination with
high sub-threshold light intensity (LI=0.75mW/mm2, blue). The slopes of the
APD90 restitution curves are 0.05, 0.06 and 0.16 for LI=0.0mW/mm2 (black),
LI=0.25mW/mm2 (red) and LI=0.75mW/mm2 (blue), respectively. The inset in (B)
shows an enlarged part of the CV restitution curve at low diastolic interval
(DI), which shows a crossover of the CV values.

uniform, global, time-periodic light pulses of 0.75mW/mm${}^{2}\leq$ LI $<$
2.0mW/mm2 led to the initiation of wave breaks, as shown in Fig. 1E. These
breaks appeared to occur right at the beginning of the ‘dark’ phase of the
applied stimulus, i.e., when the light was turned off in a pacing cycle.
Application of high subthreshold LIs led to the spontaneous emergence of
conduction blocks within the domain (see bold white lines indicated by green
arrows in Fig. 1E), which seemed to promote wave break initiation. For the
full sequence of events leading to the incidence of wave breaks, see Movie S1
of the SI Appendix.

To explain this unusual observation, we first studied the influence of light
on the dispersion of APD and CV, i.e. the restitution properties of the
system. Our results, as shown in Fig. 2A and B, pointed to a strange anomaly:
Both the APD and CV restitution curves for different LIs appeared to have
gentle slopes (<1.0) without spatial dispersion, which is fascinating because
this type of wave break intiation is typically associated with steep
restitution curves (slope > 1.0), or dispersion in restitution properties
within a system (contrary to our case), or in the presence of local ionic
heterogeneities (as in case of ischemia) (37, 38) or in domains with
progressively degrading excitability, which results in the gradual flatenning
of the CV restitution curve with reduced mean CV, characteristic of type-II
ventricular fibrillation. In addition, constant application of light to the
domain resulted in supernormal CVs, which was in direct conflict with our
observation of the emergence of lines of propagation block, just prior to the
occurrence of the break.

To understand the basis for this anomaly, we tried to reproduce the appearance
of wave blocks in a reduced system, i.e. pseudo-1D, using a rectangular
simulation domain containing 512$\times$10 grid points. We applied high
frequency electrical stimulation to the left edge of the domain (Fig. 3A),
while illuminating the entire domain uniformly, in a time-periodic manner. In
the absence of illumination, each applied stimulus resulted in the initiation
of a new wave that propagated uninterruptedly to the right edge of the domain
(Fig. 3A).

Figure 3: A) Right propagating plane waves in a pseudo-1D domain containing
512 $\times$ 10 grid points, electrically paced at 11.5 Hz from the left edge
of the domain in the absence of illumination. Each stimulus results in the
generation of a new propagating wave. B) Wave block at the pacing site when
the simulation in (A) is performed in the presence of a uniform global time-
periodic illumination, at 2.1 Hz temporal frequency and and LI=0.75 mW/mm2. C)
Wave block away from the site of pacing, in a pseudo-1D channel (1536
$\times$30 grid points) connected to a 2D domain (512 $\times$512 grid points)
that maintains a spiral. The spiral drives the electrical stimulation in the
channel. Time-periodic light stimulation is only applied to the pseudo-1D
extension of the 2D domain. Green arrows indicate the locations of the wave
block.

However, in the presence of an optical stimulation with LI$\geq$ 0.75 mW/mm2
and any frequency between 2 - 12 Hz, we observed a random blocking of the
stimulated waves, always at the site of stimulation ( Fig. 3B). Thus, based on
these data, it was impossible to conclude whether the block was actually a
dynamic effect caused by optical stimulation or simply a result of CV
restitution.

Given that wave block at the site of stimulation could give no clear insight
into the formation of zero CV lines within the domain in 2D, we shifted the
region for optical perturbation to 98 < x < 512, and resumed periodic
electrical pacing at 11.5Hz (dominant frequency of the spiral wave in the
absence of illumination), from the left end (x < 4) of the stripe. Our studies
showed that every stimulus was capable of producing a new wave that propagated
all the way to the right end of the stripe without interruption, although 1D
traces along the length of the stripe did give us signatures of the effect of
subthreshold stimulation on the membrane potential in the illuminated region.
Thus, clearly, the 1D reduction of our 2D system was not optimised for wave
block visualisation. To overcome this challenge, keeping in mind that the
spiral in our 2D system was actually a meandering one, rotating quasi-
periodically with 2 dominant frequencies, we chose to use the spiral in the 2D
domain to drive the electrical activity in the stripe domain, while applying
the optical perturbation only to the stripe domain. This had no direct bearing
on the dynamics of the spiral wave itself, but was used as a driver of the
electrical activity in the stripe domain, equivalent to the spiral in the 2D
domain. Thus, we constructed a hybrid pseudo-1D equivalent of our original 2D
system by connecting a 2048$\times$30 cell grid (pseudo-1D) to our 2D domain,
which contained 512$\times$512 grid points. We initiated a spiral in the 2D
domain. The waves emitted from the spiral core were able to move into the
pseudo-1D lattice and thereby control the frequency of the applied electrical
stimulation. We applied uniform global optical stimulation at LI=0.75mW/mm2
and low frequencies ($\simeq$ 2.0HZ) only to the pseudo-1D domain, leaving the
2D domain unaffected by light. This setup allowed us to observe the formation
of conduction blocks within the pseudo-1D domain, i.e. away from the
stimulation site, over a wide range of frequencies of the applied light
stimulation (Fig. 3C, and Movie S2 of the SI Appendix).

Figure 4: (A-G) Effect of applying a pulse of LI = 0.75 mW/mm2 and duration
100 ms to different phases of an action potential (AP). For comparison, a
control AP, i.e. in the absence of the light pulse, is provided next to the
affected AP in each sub-figure. Black curves represent the transmembrane
voltage (V), red curves show the corresponding channelrhodopsin current
($I_{ChR2}$), and yellow curves mark the duration of the application of the
light pulse.

In addition to the disappearance of some of the waves that had propagated into
the pseudo-1D domain over a certain distance, we observed spatio-temporal
oscillation in the wavelength of successive waves. To understand the basis for
these oscillations and to correlate to the appearance of the conduction block,
we studied the effects of sub-threshold illumination on a single heart cell
when exposed to light at different stages of an electrically stimulated AP
(Fig. 4). Our studies showed that for an optical stimulus that had
approximately the duration of one AP, sub-threshold illumination had the
greatest effect on the APD at $90$% repolarization (i.e. APD90), even when the
light was applied during the recovery phase of the AP (Fig. 4G).

Figure 5: (A-C) Effect of applying a 100 ms, 0.75 mW/mm2 light pulse to a
plane wave propagating through a pseudo-1D domain containing 512 $\times$ 10
points. Black and blue curves represent $V$ and $I_{ChR2}$ respectively,
measured from the line y=5 within the domain. (D) Rapid recovery of the wave
tail (indicated by bold red arrows) immediately after the light is turned off.
(E-F) Gradual recovery of the original waveform with increase in wavelength
$\lambda$. In each sub-figure, the upper (colored) panel shows the spatial
profile of the wave as it propagates through the pseudo-1D domain.

A direct consequence of this differential response of the cell membrane to
light applied at different phases of the AP is a phase-dependent delay of
recovery, that is induced along the length of a wave, in extended media. Such
a delay leads to differences in CV along a wavelength. Especially in the
presence of light, the wavefront propagates at supernormal speeds while the
waveback slows down. This results in (i) prolongation of the wavelength with a
gradual change in the shape of wave profile and (ii) blocking of the following
wave by the preceding one. When light is applied to a propagating wave, cells
located at increasing distances from the wavefront towards the waveback,
receive light at progressively shifted phases of their AP (Fig. 4A to G). This
results in a gradual decrease in the rate of recovery of the cells in between
the wavefront and the waveback, thereby changing spatial profile of the wave,
to give it a pinched appearance. (Fig. 5A-C). Instead of the constant,
uniform, global illumination, if the light is now applied periodically in
time, the pinched wavelength gets the opportunity to return to its original
state during the light-off period (Fig. 5D-F). This explains the occurrence of
spatiotemporal oscillations in the wavelength of the propagating wave (in 1D)
or a spiral arm (in 2D).

Figure 6: Demonstration of the mechanism of spiral wave breakup using a series
of 1D traces of the spatial profile of the electrical activity in the medium,
as measured along the solid white arrow shown in the inset. Here, y is the
displacement from the center of the domain to the boundary, along the length
of the arrow. The time stamp for each 1D trace is indicated on the left hand
side of the vertical axis. The domain is illuminated uniformly and globally
with LI = 0.75 mW/mm2, from 430 ms to 610 ms. During this time, the CV of the
right-propagating wavefront increases beyond normal, whereas, that of the
waveback decreases. This manifests in a change in the spatial wave profile
during the said period. Towards 570 ms, the wave profile modifies sufficiently
to allow creation of a vulnerable window. Subsequent removal of light at 610
ms opens the window for re-excitation, which then permits re-entry; re-
excitation occurs at 620 ms, which then gives rise to a left-propagating wave.
The translucent pink rectangle marks the spatiotemporal window for wave break
initiation. Reversal of wave direction is indicated on the figure using a
colour gradient from grey, for left-propagating waves, to white, for
propagation to the right.

The higher the LI, the stronger the conditioning of the wavelength, i.e. the
wave appears more pinched. The moment the light is switched off, the
differences between the recovery rates of the cells at different distances
from the wavefront disappear. The cells on the wavefront that allowed
propagation at supernormal CV recover slowly, allowing the wavefront to expand
from its pinched state in a gradual manner. The cells on the wave back that
were forcibly held at a sub-threshold voltage higher than the normal resting
membrane potential (RMP) in the pinched state recover immediately, creating a
window for re-excitation. In 2D, this window allows re-entry, as shown in Fig.
6. The inset shows the spatial pattern of electrical activity in the domain,
which assumes a spiral form. By tracing the voltage distribution along the
solid white line (see Fig. 6, inset) in the direction of the arrow, we
illustrate how a right-propagating wave is modulated over time with an optical
perturbation (LI = 0.75 mW/mm2, from 430 ms to 610 ms) to create a vulnerable
window (marked with a pink rectangle). The wave break occurs between 600ms and
620ms, leading to reexcitation and reentry, marked by the initiation of a
left-propagating wave. Reversal of wave direction is indicated on the figure
using a colour gradient from grey, for left-propagating waves, to white, for
propagation to the right.

Finally, to test the generic nature of this new mechanism of wave break
initiation at stable parameter régimes, we repeated the study in another
system, a 2D model of neonatal mouse ventricular tissue. At uniform global
constant illumination with LI = 0.025 mW/mm2, a spiral wave rotated with no
signatures of breakup (Fig. 7A). However, when stimulated periodically at 4
Hz, with light of the same intensity and pulse length 200 ms, breakup similar
to Fig. 1E was observed (Fig. 7B). A study of the restitution characteristics
for this system, as presented in Ref.(39) (see Figures 5 A and B, therein),
indicated that wave breaks occurred in a stable parameter régime, thus proving
the robustness and restitution-independence of our proposed mechanism for wave
breaks.

Figure 7: Spiral wave breakup in the 2D model for optogenetically-modified
neonatal mouse ventricular tissue. A) Spiral wave in the presence of uniform,
global, constant illumination at LI = 0.025mW/mm2. B) Break up of spiral waves
in the presence of periodic optical stimulation at 4 Hz, with a light pulse of
duration 200 ms.

##  2 Discussion

An intriguing question is: Would the choice of model parameters have an
influence on wave break initiation, through a parameter-induced destabilizing
effect (which may be characteristic of the cAF model used)? The answer is no.
As such, the cAF model plays no specific role in determining the mechanism of
wave break initiation and reentry. In Fig. 7 we demonstrate using a neonatal
mouse ventricular model that such breaks can also occur in the ventricle. The
reason for using the cAF model is because (i) it makes more sense to
‘defibrillate’ a diseased condition rather than the healthy heart, where
fibrillation is less likely to occur. We believe that any other parameter
régime that supports fibrillation (atrial or ventricular type) or tachycardia,
should also demonstrate this phenomenon, under the application of pulsed
subthreshold light stimulation. (ii) The human atrial wall is super thin (
2-3mm) in contrast to the ventricles ( 1.5cm), which makes the former a better
candidate for optogenetic defibrillation, whose efficiency seems to rely on
maximal illumination of bulk of the affected tissue. Researchers believe that
the problem of poor penetration of visible light in cardiac tissue might not
be so important in the human atria, because of its remarkably low wall
thickness. Nevertheless, our studies reveal that one needs to be careful while
designing the defibrillation protocol, as a series of globally illuminating
light pulses in the suprathreshold regime, can, in fact trigger the formation
of wave breaks underneath the surface, where illumination is subthreshold, and
thereby lead to failure of defibrillation. Finally, (iii) from a technical
perspective, we chose the cAF model over the healthy heart model because of
the large computational costs associated with the latter (see Methods for
details).

There are several theories regarding the factors and conditions that cause a
spiral wave to break up. Of these, the strongest candidate is restitution
(both APD and CV) in cardiac tissue (19). Some studies argue that APD and CV
restitution curves fundamentally characterize the wave dynamics in the heart
as they reflect the mesoscopic effects of changes in ion currents and
concentrations occurring at the cellular level (19). Typically, a steep APD
restitution curve (slope $\geq 1.0$) can drive the system via Hopf bifurcation
to APD oscillations (40, 41, 42). Such alternans promotes the development of
functional conduction blocks because the propagation of a "long" action
potential (AP) into a "short" AP region fails, thereby limiting the allowable
diastolic interval (DI) (19). Another proposed mechanism, again relating to
the role of a steep APD restitution curve, suggests that the steepness of the
curve is negatively correlated with the speed of the waveback in the presence
of recovery gradients (43). A gradient $\geq 1.0$ introduces a difference
between the CVs of the propagating wavefront and back so that within a wave
train the front of the following wave collides with the back of the preceding
wave, resulting in a functional conduction block (41, 44, 43, 45). Other
mechanisms of wave break include the occurrence of a spatially non-
synchronized type of APD oscillation, also known as spatially discordant APD
alternans (46, 47, 48, 49), the presence of a biphasic APD restitution curve
(50, 51), effects of hysteresis (52, 53, 54, 55) and Doppler shift through the
trajectories of the spiral peaks (56).

A fascinating study by Fenton et al. (57) on wave break initiation in ischemic
heart tissue shows that dynamic flatenning of the APD restitution curve
(characterized by slope < 1.0), results in an electrical ‘memory’ effect,
which provokes a monotonic reduction in the excitability of cardiac tissue.
Over time, the excitability becomes so low that a spiral breaks up close to
its tip in the form of a ‘voltage drop’. The study concludes that spiral
break-up can be induced in parameter régimes with flat APD restitution as long
as the dynamic effects of tissue memory and substantially large changes in
basic cycle length are taken into account. Our study is consistent with the
work of Fenton et al. in that wave breaks are observed at low APD and CV
restitution slopes. However, we show that for a system with static APD and CV
restitution curves, in the absence of a factor leading to a progressive
reduction in tissue excitability and without explicit inclusion of the
electrical memory effect, it is still possible to initiate wave breaks that
develop into sustained re-entry.

Another study by Zemlin et al. (58) reports spiral wave break-up in parameter
regimes characterized by strongly negative slopes of the APD restitution
curve. However, as noted in their work, such break-up is often less pronounced
than the break-up that occurs with positive slopes, mainly because break-up
with negative restitution takes a fairly long time to mature, compared to
break-up with positive restitution. The wave break initiation that we
demonstrate in this article happens very fast, in most cases, requiring not
more than 1-2 pulses.

Another study worth mentioning in this context is that of Banville and Gray
(59), which demonstrates in experiments on the rabbit heart, the crucial role
of the spatial distribution of APD- and CV restitution on the onset of
alternans and fibrillation. However, in our study, we used a homogeneous
domain with all cells identical, and subjected each cell to the same external
conditions. Thus we did not have a dispersion of APD and CV restitution that
could influence wave break initiation like in the case of (59). Nevertheless,
we observe the initiation of wave breaks, which leads us to rethink the
underlying factors involved in initiating wave breaks in cardiac tissue under
different conditions.

A rare, but quite efficient, mechanism of spiral wave break up involves
supernormal excitability in certain domains. It is associated with a decrease
in the excitation threshold with decreasing diastolic interval (DI).
Typically, if it exists in a model, it promotes the development of conduction
blocks through collisions between the rapidly advancing wave fronts and
subsequent stacking at short DI, eventually leading to the initiation of wave
breaks. Breakup occurs particularly close to the spiral tip. Supernormal CV
can also lead to heterogeneity of refractoriness at shallow APD restitution,
in which case scalloping occurs. The mechanism of spiral wave break initiation
that we report in the present article is somewhat similar to this phenomenon.
With uniform global light perturbations, we are able to impose supernormal CVs
on the propagating electrical activity. This results in the development of
wave blocks within the domain that are visible as white lines in Fig. 1E. In
our studies, the reentry occurs once the factor causing supernormality is
withdrawn, i.e., the light perturbation is removed, allowing the wavefront and
wave back to restore their normal velocities.

Our study provides a clear example of a funny mechanism of wave break
initiation and reentry that occurs when the slope of the APD restitution curve
is less than 1.0 and the CV restitution curve is relatively flat. The proposed
mechanism should be looked on as an important factor that can potentially lead
to failure of optical defibrillation using pulsed global illumination. We
believe that this mechanism is important to note because it provides useful
insights into the translational prospects of optogenetic defibrillation in the
human atria. And as our study proves, the phenomenon is non crucially
dependent on the choice of the model used. Thus, we provides new insights into
the possibility of the occurrence of wave breaks in heart tissue, even when
there are no standard markers for vulnerability.

Finally, the effect of global subthreshold stimulation was studied in simple
ionic models of cardiac tissue, in the context of electrical defibrillation. A
particularly close phenomenon was reported by Sridhar and Sinha (60), who
demonstrated the possibility to terminate electrical activity using a time-
invariant global subthreshold electrical stimulation to achieve
synchronization. Their study indicated that an electric pulse of length
$\mathcal{O}(1APD)$ was sufficient to cause this effect. Our studies showed
(data not presented) that if the stimulus in(60) was instead, applied at low
frequency ( 1-2Hz), using pulses of very short duration (
$\mathcal{O}(0.1APD)$), then such perturbation could also initiate wavebreaks.
However, such breaks occur within a very restricted parameter régime and
critically rely on the frequency of the applied perturbation for sustenance,
owing to the sensitive dependence of the response of the single cell AP
morphology, on the phase of the AP at which the stimulus is applied.

Electrical activity in cardiac tissue was modeled using the following
reaction-diffusion-type equation:

$\frac{dV}{dt}=\nabla.\mathcal{D}\nabla V-\frac{I_{ion}}{C_{m}}$ (1)

where $V$ represents the transmembrane potential (in mV) developed across
single cardiomyocytes, $C_{m}$ is the specific capacitance (in ${\rm\mu
F/cm^{2}}$) of the cell membrane, $\mathcal{D}$ represents the diffusion
coefficient for intercellular coupling, and $I_{ion}$ represents the total
ionic current produced by each individual cell. For human atrial tissue,
$I_{ion}$ was formulated according to the Courtemanche-Ramirez-Nattel (CRN)
model (61). It is a sum of $12$ ionic currents: fast $Na^{+}$ ($I_{Na}$),
inward rectifier $K^{+}$ ($I_{K1}$), transient outward $K^{+}$ ($I_{to}$),
ultra-rapid $K^{+}$ ($I_{Kur}$), rapid and slow delayed rectifier $K^{+}$
($I_{Kr}$ and $I_{Ks}$, respectively), $Na^{+}$ and $Ca^{2+}$ background
($I_{BNa}$, and $I_{CaL}$), $Na^{+}/K^{+}$ pump, $Ca^{2+}$ pump ($I_{pCa}$),
$Na^{+}/Ca^{2+}$ exchanger ($I_{NaCa}$) and $L-type\leavevmode\nobreak\
Ca^{2+}$ current ($I_{CaL}$). The parameters of the CRN model were adjusted to
reproduce the action potential of the normal atrial working myocardium (62).
In particular, the maximum conductance for $I_{K1}$, i.e., $G_{K1}$ was
changed from 0.09 nS/pF in the original model, to 0.117 nS/pF, to yield a
single cell action potential duration at 90% repolarization of the membrane
potential (APD90), of 284 ms. A choice of $\mathcal{D}$=0.0023 cm2/ms,
produced a conduction velocity (CV) of 69.75 cm/s in the two-dimensional
healthy tissue domain. In order to model the action potential during chronic
atrial fibrillation (AF) remodelling, the maximal conductances of $I_{to}$,
and $I_{CaL}$ were reduced by $85\%$, and $74\%$, respectively, $G_{K1}$ was
increased by $250\%$, the time constant for activation of $I_{CaL}$ was
increased by $62\%$, the activation curves for $I_{to}$ and $I_{CaL}$ wer
shifted by +16 mV and -5.4 mV, respectively, while the inactivation curve for
$I_{Na}$ was shifted by +1.6 mV (63, 64, 65). These parameter adjustments
reduced the wavelength (at 1Hz electrical pacing) from $\simeq$18 cm in
healthy tissue, to $\simeq$5 cm, allowing us to fit the spiral into a smaller
simulation domain (512 $\times$ 512, as opposed to 2048 $\times$ 2048) without
causing it to break up. This improved the computational cost-effectiveness by
a factor of 20. The spiral wave in our study meandered with a hypocycloidal
tip trajectory, and survived for longer than 10s of simulation time.

For the neonatal mouse ventricular model, $I_{ion}$ was expressed as a sum of
16 ionic currents, according to the Wang and Sobie (66). These include, the
fast $Na^{+}$ current ($I_{Na}$), the background $Na^{+}$ and $Ca^{2+}$
currents ($I_{Nab}$ and $I_{Cab}$), the L-type and T-type $Ca^{2+}$ currents
($I_{CaL}$ and $I_{CaT}$), the $Ca^{2+}$ pump current ($I_{pCa}$), the
$Na^{+}/Ca^{2+}$ exchanger ($I_{NaCa}$), $Na^{+}/K^{+}$ pump ($I_{NaK}$), slow
and fast components of the transient outward $K^{+}$ currents ($I_{Kto,f}$ and
$I_{Kto,s}$), slow and rapid delayed rectifier $K^{+}$ currents ($I_{Ks}$ and
$I_{Kr}$), ultrarapid delayed rectifier $K^{+}$ current ($I_{Kur}$), sustained
outward $K^{+}$ current ($I_{Kss}$), inward rectifier ($I_{K1}$), and the
$Ca^{2+}$-activated $Cl^{-}$ current ($I_{Cl,Ca}$). We used
$\mathcal{D}$=0.00095 cm/ms, which led to a CV of 43.9 cm/s. The 2D
simnulation domain contained 200 $\times$ 200 grid points.

In order to incorporate the effects of optogenetics, our models were combined
with a 4-state model for voltage- and light-sensitive Channelrhodopsin-2 (67),
as described in Eqs. 2-7.

$\displaystyle I_{ChR2}=g_{ChR2}G(V)(O_{1}+\gamma O_{2})(V-E_{ChR2},$ (2)
$\displaystyle\frac{dC_{1}}{dt}=G_{r}C_{2}+G_{d1}O_{1}-k_{1}C_{1},$ (3)
$\displaystyle\frac{dC_{2}}{dt}=G_{d2}O_{2}-(k_{2}+G_{r})C_{2},$ (4)
$\displaystyle\frac{dO_{1}}{dt}=k_{1}C_{1}-(G_{d1}+e_{12})O_{1}+e_{21}O_{2},$
(5)
$\displaystyle\frac{dO_{2}}{dt}=k_{2}C_{2}-(G_{d2}+e_{21})O_{2}+e_{12}O_{1},$
(6) $\displaystyle O_{1}+O_{2}+C_{1}+C_{2}=1$ (7)

Here $O_{1}$, $O_{2}$, $C_{1}$ and $C_{2}$ represent the open and closed
states of the ChR2 ion channel. $k_{1}$, $k_{2}$, $G_{d1}$, $G_{d2}$, $G_{r}$,
$e_{12}$, and $e_{21}$ are the kinetic parameters, whereas, $G_{r}$, $G(V)$
and $G_{d1}$ contain the voltage-dependencies. For a detailed description of
these parameters and their values, we refer the reader to (67). The total
current produced by each cell in Eq. 1 is $I_{ion}$, plus, the
Channelrhodopsin current $I_{ChR2}$. We found that for human atrial tissue, LI
< 2.0mW/mm2 were typically sub-threshold, meaning, they failed to trigger
action potentials, whereas, higher light intensities successfully stimulated
cells and initiated waves in extended media. In neonatal mouse, LI < 0.03
mW/mm2 was considered sub-threshold.

This work was supported by the Max Planck Society and the German Center for
Cardiovascular Research (DZHK).

## SI Appendix

Movie S1: Controlled break-up of spiral waves, induced by periodic sub-
threshold optical stimulation, at frequency 1.56 Hz. The trajectory of the
spiral tip(s) is traced using white lines that fade over time. The filled blue
circle below the voltage color-code (in mV) in select frames indicates the
application of uniform, global subthreshold illumination at LI = 0.75 mW/mm2,
in these frames.

Movie S2: Reduction of the 2D system to a hybrid 1D system, in which we
combine a 2D simulation domain to a pseudo-1D outlet. The electrical
stimulation in the pseudo-1D domain is driven by the unperturbed spiral wave
in 2D. A uniform global periodic sub-threshold stimulation (frequency: 2.1 Hz)
with light of intensity 0.75 mW/mm2 is applied to the pseudo-1D domain.
Propagation failure occurs well within the pseudo-1D domain for randomly
selected waves. In addition, we observe oscillations in wavelength and
modulation of the wave profile for waves propagating through the pseudo-1D
domain.

## References

  * (1) E Pálsson, EC Cox, Origin and evolution of circular waves and spirals in dictyostelium discoideum territories.  Proc Natl Acad Sci U S A 93, 1151–1155 (1996).
  * (2) RJ Field, M Burger, eds., Oscillations and Traveling Waves in Chemical Systems. (Wiley, New York), (1985).
  * (3) AM Zhabotinsky, A history of chemical oscillations and waves.  Chaos 1, 379–86 (1991).
  * (4) M Bär, N Gottschalk, M Eiswirth, G Ertl, Spiral waves in a surface reaction: Model calculations.  The Journal of Chemical Physics 100, 1202–1214 (1994).
  * (5) V Krinsky, H Swinney, eds., Wave and Patterns in Biological and Chemical Excitable Media.  (North-Holland, Amsterdam), (1991).
  * (6) R Kapral, K Showalter, eds., Chemical Waves and Patterns.  (Kluwer, Dordrecht), (1995).
  * (7) JM Davidenko, P Kent, DR Chialvo, DC Michaels, J Jalife, Sustained vortex-like waves in normal isolated ventricular muscle.  Proc. Natl. Acad. Sci. USA. 87, 8785–8789 (1990).
  * (8) AV Panfilov, RH Keldermann, MP Nash, Drift and breakup of spiral waves in reaction–diffusion–mechanics systems.  Proceedings of the National Academy of Sciences 104, 7922–7926 (2007).
  * (9) NF Otani, A primary mechanism for spiral wave meandering.  Chaos: An Interdisciplinary Journal of Nonlinear Science 12, 829–842 (2002).
  * (10) H Zhang, A Holden, Chaotic meander of spiral waves in the fitzhugh-nagumo system.  Chaos, Solitons & Fractals 5, 661 – 670 (1995) Nonlinear Phenomena in Excitable Physiological Systems.
  * (11) O Steinbock, SC Müller, Light-controlled anchoring of meandering spiral waves.  Phys. Rev. E 47, 1506–1509 (1993).
  * (12) P Kuklik, JJ Zebrowski, L Szumowski, P Sanders, Anchoring of spiral waves by inhomogeneity of anisotropic conduction.  Heart, Lung and Circulation 17, S7–S8 (2008).
  * (13) P Boily, VG LeBlanc, E Matsui, Spiral anchoring in media with multiple inhomogeneities: A dynamical system approach.  Journal of Nonlinear Science 17, 399–427 (2007).
  * (14) D Olmos, Reflection and attachment of spirals at obstacles for the fitzhugh-nagumo and beeler-reuter models.  Phys. Rev. E 81, 041924 (2010).
  * (15) H Kitahata, M Tanaka, Mathematical approach to unpinning of spiral waves anchored to an obstacle with high-frequency pacing.  Biophysics and physicobiology 15, 196–203 (2018).
  * (16) P Porjai, et al., Electrically forced unpinning of spiral waves from circular and rectangular obstacles.  Chemical Physics Letters 660, 283 – 286 (2016).
  * (17) X Feng, X Gao, DB Pan, BW Li, H Zhang, Unpinning of rotating spiral waves in cardiac tissues by circularly polarized electric fields.  Scientific Reports 4 (2014).
  * (18) M Tanaka, et al., Unpinning of a spiral wave anchored around a circular obstacle by an external wave train: Common aspects of a chemical reaction and cardiomyocyte tissue.  Chaos: An Interdisciplinary Journal of Nonlinear Science 19, 043114 (2009).
  * (19) FH Fenton, EM Cherry, HM Hastings, SJ Evans, Multiple mechanisms of spiral wave breakup in a model of cardiac electrical activity.  Chaos: An Interdisciplinary Journal of Nonlinear Science 12, 852–892 (2002).
  * (20) S Luther, et al., Low-energy control of electrical turbulence in the heart.  Nature 475, 235 – 239 (2011).
  * (21) L Xu, Z Li, Z Qu, Z Di, Resonance drifts of spiral waves on media of periodic excitability.  Physical review. E, Statistical, nonlinear, and soft matter physics 85, 046216 (2012).
  * (22) S Grill, VS Zykov, SC Müller, Spiral wave dynamics under pulsatory modulation of excitability.  The Journal of Physical Chemistry 100, 19082 – 19088 (1996).
  * (23) J Ramos, Dynamics of spiral waves in excitable media with local time-periodic modulation.  Chaos, Solitons & Fractals 13, 1383 – 1392 (2002).
  * (24) O Steinbock, V Zykov, S Müller, Control of spiral-wave dynamics in active media by periodic modulation of excitability.  Nature 366, 322–324 (1993).
  * (25) A Schrader, M Braune, H Engel, Dynamics of spiral waves in excitable media subjected to external periodic forcing.  Phys. Rev. E 52, 98–108 (1995).
  * (26) J Ma, CN Wang, J Tang, YF Xia, Suppression of the spiral wave and turbulence in the excitability-modulated media.  International Journal of Theoretical Physics 48, 150 – 157 (2009).
  * (27) T Bruegmann, et al., Optogenetic control of heart muscle in vitro and in vivo.  Nat Methods. 7, 897–900 (2010).
  * (28) C Crocini, et al., Optogenetics design of mechanistically-based stimulation patterns for cardiac defibrillation.  Sci Rep. 6, 35628 (2016).
  * (29) T Bruegmann, et al., Optogenetic defibrillation terminates ventricular arrhythmia in mouse hearts and human simulations.  The Journal of Clinical Investigation 126, 3894–3904 (2016).
  * (30) ECA Nyns, et al., Optogenetic termination of ventricular arrhythmias in the whole heart: towards biological cardiac rhythm management.  Eur Heart J. 38, 2132–2136 (2017).
  * (31) T Bruegmann, T Beiert, C Vogt, J Schrickel, P Sasse, Optogenetic termination of atrial fibrillation in mice.  Cardiovasc Res. 114, 713–723 (2018).
  * (32) R Majumder, et al., Optogenetics enables real-time spatiotemporal control over spiral wave dynamics in an excitable cardiac system.  eLife 7, e41076 (2018).
  * (33) Z Jia, et al., Stimulating cardiac muscle by light: cardiac optogenetics by cell delivery.  Circ Arrhythm Electrophysiol. 4, 753–60 (2011).
  * (34) I Feola, et al., Localized optogenetic targeting of rotors in atrial cardiomyocyte monolayers.  Circ Arrhythm Electrophysiol. 10, e005591 (2017).
  * (35) S Park, S Lee, L Tung, D Yue, Optical mapping of optogenetically shaped cardiac action potentials.  Scientific reports 4, 6125 (2014).
  * (36) R Burton, et al., Optical control of excitation waves in cardiac tissue.  Nature photonics 9, 813 (year?).
  * (37) T Wu, S Lin, J Weiss, C Ting, P Chen.  Circulation 106, 1859–1866 (2002).
  * (38) Y Liu, et al.  J. Cardiovasc. Electrophysiol. 15, 1433–1440 (2004).
  * (39) A Mayer, P Bittihn, S Luther, Complex restitution behavior and reentry in a cardiac tissue model for neonatal mice.  Physiological Reports 5, e13449 (2017).
  * (40) JB Nolasco, RW Dahlen, A graphic method for the study of alternation in cardiac action potentials.  Journal of Applied Physiology 25, 191–196 (1968).
  * (41) M Courtemanche, L Glass, JP Keener, Instabilities of a propagating pulse in a ring of excitable media.  Phys. Rev. Lett. 70, 2182–2185 (1993).
  * (42) A Karma, H Levine, X Zou, Theory of pulse instabilities in electrophysiological models of excitable tissues.  Physica D: Nonlinear Phenomena 73, 113 – 127 (1994).
  * (43) M Courtemanche, Complex spiral wave dynamics in a spatially distributed ionic model of cardiac electrical activity.  Chaos: An Interdisciplinary Journal of Nonlinear Science 6, 579–600 (1996).
  * (44) M Courtemanche, JP Keener, L Glass, A delay equation representation of pulse circulation on a ring in excitable media.  SIAM J. Appl. Math. 56, 119–142 (1996).
  * (45) M Watanabe, F Fenton, S Evans, H Hastings, A Karma, Mechanisms for discordant alternans.  Journal of Cardiovascular Electrophysiology 12, 196–206 (2001).
  * (46) DS Rosenbaum, et al., Electrical alternans and vulnerability to ventricular arrhythmias.  New England Journal of Medicine 330, 235–241 (1994).
  * (47) D Sato, et al., Spatially discordant alternans in cardiac tissue.  Circulation Research 99, 520–527 (2006).
  * (48) H Hayashi, et al., Dynamic origin of spatially discordant alternans in cardiac tissue.  Biophysical journal 92, 448–460 (2007).
  * (49) R Majumder, MC Engels, AF de Vries, Antoine, AV Panfilov, A Pijnappels, Daniel, Islands of spatially discordant apd alternans underlie arrhythmogenesis by promoting electrotonic dyssynchrony in models of fibrotic rat ventricular myocardium.  Scientific reports 6 (2016).
  * (50) MR Franz, CD Swerdlow, LB Liem, J Schaefer, Cycle length dependence of human action potential duration in vivo. effects of single extrastimuli, sudden sustained rate acceleration and deceleration, and different steady-state frequencies.  The Journal of Clinical Investigation 82, 972–979 (1988).
  * (51) P Szigligeti, et al., Intracellular calcium and electrical restitution in mammalian cardiac cells.  Acta Physiologica Scandinavica 163, 139–147 (1998).
  * (52) P Lorente, J Davidenko, Hysteresis phenomena in excitable cardiac tissuesa.  Annals of the New York Academy of Sciences 591, 109–127 (1990).
  * (53) GM Hall, S Bahar, DJ Gauthier, Prevalence of Rate-Dependent Behaviors in Cardiac Muscle.  Physical Review Letters 82, 2995–2998 (1999).
  * (54) AR Yehia, D Jeandupeux, F Alonso, MR Guevara, Hysteresis and bistability in the direct transition from 1:1 to 2:1 rhythm in periodically driven single ventricular cells.  Chaos: An Interdisciplinary Journal of Nonlinear Science 9, 916–931 (1999).
  * (55) RA Oliver, et al., Existence of bistability and correlation with arrhythmogenesis in paced sheep atria.  Journal of Cardiovascular Electrophysiology 11, 797–805 (2000).
  * (56) O Steinbock, J Schütze, SC Müller, Electric-field-induced drift and deformation of spiral waves in an excitable medium.  Phys. Rev. Lett. 68, 248–251 (1992).
  * (57) F Fenton, S Evans, H Hastings, Memory in an excitable medium: a mechanism for spiral wave breakup in the low-excitability limit.  Physical review letters 83, 3964 (1999).
  * (58) C Zemlin, A Panfilov, Spiral waves in excitable media with negative restitution.  Physical Review E. 63, 041912 (2001).
  * (59) I Banville, R Gray, Effect of action potential duration and conduction velocity restitution and their spatial dispersion on alternans and the stability of arrhythmias.  Journal of cardiovascular electrophysiology 13, 1141–9 (2002).
  * (60) S Sridhar, S Sinha.  Euro Physics Letters 92, 60006 (2010).
  * (61) M Courtemanche, RJ Ramirez, S Nattel, Ionic mechanisms underlying human atrial action potential properties: insights from a mathematical model.  American Journal of Physiology-Heart and Circulatory Physiology 275, H301–H321 (1998) PMID: 29586616.
  * (62) C Tobón, et al., A three-dimensional human atrial model with fiber orientation. electrograms and arrhythmic activation patterns relationship.  PLOS ONE 8, 1–13 (2013).
  * (63) M Courtemanche, RJ Ramirez, S Nattel, Ionic targets for drug therapy and atrial fibrillation-induced electrical remodeling: insights from a mathematical model.  Cardiovascular Research 42, 477–489 (1999).
  * (64) H Zhang, CJ Garratt, J Zhu, AV Holden, Role of up-regulation of IK1 in action potential shortening associated with atrial fibrillation in humans.  Cardiovascular Research 66, 493–502 (2005).
  * (65) SV Pandit, et al., Ionic determinants of functional reentry in a 2-d model of human atrial cells during simulated chronic atrial fibrillation.  Biophysical Journal 88, 3806 – 3821 (2005).
  * (66) LJ Wang, EA Sobie, Mathematical model of the neonatal mouse ventricular action potential.  American Journal of Physiology-Heart and Circulatory Physiology 294, H2565–H2575 (2008) PMID: 18408122.
  * (67) JC Williams, et al., Computational optogenetics: Empirically-derived voltage- and light-sensitive channelrhodopsin-2 model.  PLOS Computational Biology 9, 1–19 (2013).


# Area Efficient Modular Reduction in Hardware for Arbitrary Static Moduli

Robin Müller, Willi Meier, Christoph F. Wildfeuer
FHNW University of Applied Sciences and Arts Northwestern Switzerland
School of Engineering
5210 Windisch, Switzerland
<EMAIL_ADDRESS>
###### Abstract

Modular reduction is a crucial operation in many post-quantum cryptographic
schemes [1, 2], including the Kyber key exchange method [3] or Dilithium
signature scheme [4]. However, it can be computationally expensive and pose a
performance bottleneck in hardware implementations. To address this issue, we
propose a novel approach for computing modular reduction efficiently in
hardware for arbitrary static moduli. Unlike other commonly used methods such
as Barrett or Montgomery reduction, the method does not require any
multiplications. It is not dependent on properties of any particular choice of
modulus for good performance and low area consumption. Its major strength lies
in its low area consumption, which was reduced by 60% for optimized and up to
90% for generic Barrett implementations for Kyber and Dilithium. Additionally,
it is well suited for parallelization and pipelining and scales linearly in
hardware resource consumption with increasing operation width. All operations
can be performed in the bit-width of the modulus, rather than the size of the
number being reduced. This shortens carry chains and allows for faster
clocking. Moreover, our method can be executed in constant time, which is
essential for cryptography applications where timing attacks can be used to
obtain information about the secret key.

###### Index Terms:

Modular Reduction, FPGA, Digital Hardware, Cryptography

##  I Introduction

Computation of the modulo operation is an important part of many cryptographic
schemes, where operations are performed in a finite field. Among these are
multiple candidates for post-quantum cryptographic schemes [1], such as Kyber
(key encapsulation) [3] or Dilithium (signature) [4]. Kyber is an IND-
CCA2-secure key encapsulation mechanism, whose security is based on the
hardness of solving the learning-with-errors problem over module lattices. The
computation of the mod operation is among the most often performed operations
in these schemes. It has to be performed many times after other arithmetic
operations, such as additions or multiplications. Computing modular reductions
efficiently can therefore help improve the overall performance and efficiency
of these cryptographic schemes. There are multiple common techniques for
computing the reduction in hardware, such as Barrett reduction [5] or
Montgomery reduction [6]. A major advantage of these schemes, for software
implementations, is that they require only few steps in sequential execution.
In the case of hardware implementations, these schemes are still relatively
expensive, since they rely heavily on multiplications. These consume
significantly more chip area compared to other operations, such as additions.
A commonly used technique to make these methods more hardware friendly is to
have multiplication constants that contain relatively few ones in binary
representation, making multiplication with them feasible by a few shifts and
add operations. For this technique to be feasible, one needs to rely on the
underlying scheme to use moduli for which Barrett constants can be constructed
that contain only few binary ones. This makes the method susceptible to
different choices of the modulus. In this paper, we show a method of
performing modular reduction for a general modulus without using
multiplication, using only a mixture of additions, look-up-tables and
multiplexers. The implementation shows similar or smaller hardware use and
performance on FPGAs as a number optimized Barrett approach, but without
relying on any specific choice of the modulus.

The underlying idea of using various forms of precomputed values for reduction
has been considered before [7, 8, 9]. Most often, these attempts focus on
efficient software implementation, which can differ significantly from
efficient hardware implementations. For hardware implementations, a similar
approach to ours has been shown in [10], where the underlying idea was used to
form an iterative scheme. In this article, we present an optimized
parallelized approach to modular reduction using precomputed values. It
reduces area consumption compared to the Barrett method and increases data
throughput compared to the related iterative approach [10].

In Section II we first introduce the commonly used Barrett reduction which we
use as a reference to compare our scheme. We then introduce our method in III
and discuss the underlying idea as well as the intended hardware
implementation architecture. Finally, we show the implementation results on an
FPGA platform and perform ASIC synthesis. We compare results for both Barrett
and our approach in IV.

##  II Barrett Reduction

In this article, we compare our scheme to the commonly used Barrett reduction
[5] method. We therefore briefly introduce the Barrett reduction as well as a
commonly made optimization to it for hardware implementations. The main scheme
discussed in this article (stating from III), does not rely on the Barrett
reduction.

###  II-A General Barrett Reduction

In this description of the Barrett reduction, we let $\tilde{c}$ be a number
that we want to reduce $\mod q$ to obtain the reduced number $c$, where
$c\leftarrow\tilde{c}\mod q$. We do this by subtracting the correct multiple
of the modulus $q$ from $\tilde{c}$ to finally receive our result $c$ where
$c<q$. To find the correct multiple of $q$, we divide by a power of $2$, which
can simply be performed in hardware using shift operations. We store a
constant $m$ which approximates the ratio between the performed power of two
division and the modulus $q$. By multiplying the division result by $m$, we
obtain the (almost) correct multiple of $q$ that needs to be subtracted. In
Alg. 1 we show Barrett reduction in the general case, which was taken from
[11].

Algorithm 1 General Barrett Reduction mod $q$ [11]

1:$\tilde{c}\in[0,q^{2})$, $m$ and $k$ such that $m=\lfloor 2^{k}/q\rfloor$

2:$c=\tilde{c}\mod q$

3:$t\leftarrow(\tilde{c}\cdot m)\gg k$

4:$c\leftarrow\tilde{c}-(t\cdot q)$

5:if $c\geq q$ then

6: $c\leftarrow c-q$

7:end if

8:return $c$

###  II-B Modulus Optimized Barrett Reduction

To optimize the Barrett reduction for better hardware resource efficiency, we
mainly need to look at the multiplications in the scheme. Let $a,b$ be two
factors that we want to multiply. We let $N$ be the index (starting from $0$)
of the most significant bit in the larger of the two factors ($a$ or $b$) as
defined in (1). We can compute the multiplication of $a$ times $b$ in binary
representation, as shown in (2). For this, we shift $a$ by all integer values
starting from $0$ to $N$ and sum those terms, where $b$ contains a binary one
at the index by which we shifted $a$ (binary multiplication by either $1$ or
$0$ in (2)):

$N+1=\lceil\log_{2}{\max(\\{a,b\\})}\rceil$ (1) $a\cdot
b=\sum_{i=0}^{N}(b_{i}\cdot(a_{i}\ll i))$ (2)

If in (2) the binary representation of $b$ only contains a few ones, the
multiplication with $b$ becomes much simpler, since only a few shifted terms
of $a$ need to be considered in the sum. All terms where $b_{i}$ is zero can
be omitted. We can utilize this property to build efficient multipliers in
cases where $a$ is a variable term but $b$ is constant. We now want to choose
$b$ such that it contains as few binary ones as possible. In the case of
Barrett reduction as introduced in Alg. 1, those constant terms by which we
need to multiply, would be $q$ as well as $m$. Alg. 2 shows how this scheme
would be applied to Barrett reduction in the case of the Kyber Post Quantum
Cryptography algorithm [12], where $q=3329$. In this case, the multiplication
by $q$ benefits heavily from this optimization, since 3329 (binary:
$110100000001$) only contains 4 binary ones. Multiplication by $m$ only
benefits from it to a limited extent, since it contains 9 binary ones. While
this form of optimization can lead to highly efficient modular reduction in
some cases of $q$ and $m$, the gained resource efficiency through
optimization, can vary significantly between different choices of $q$ and $m$.
Furthermore, the value of $q$ (and $m$ indirectly since it depends on $q$) is
most often given by the cryptographic scheme and can often not be chosen
solely on the criteria of efficient Barrett reduction.

Algorithm 2 Barrett mod 3329

1:$q=2^{11}+2^{10}+2^{8}+1,m=5039,k=24,\tilde{c}\in[0,q^{2})$

2:$c=\tilde{c}\mod q$

3:$t\leftarrow(\tilde{c}\ll 12)+(\tilde{c}\ll 9)+(\tilde{c}\ll
8)+(\tilde{c}\ll 7)+(\tilde{c}\ll 5)+$

4:$(\tilde{c}\ll 3)+(\tilde{c}\ll 2)+(\tilde{c}\ll 1)+\tilde{c}$
$\triangleright$ Equiv. $m\cdot\tilde{c}$

5:$t\leftarrow t\gg k$

6:$t\leftarrow(t\ll 11)+(t\ll 10)+(t\ll 8)+t$ $\triangleright$ Equiv. $q\cdot
t$

7:$c\leftarrow\tilde{c}-t$

8:if $c\geq q$ then

9: $c\leftarrow c-q$

10:end if

11:return $c$

##  III Algorithm

In this section, we show the proposed method for computing modular reduction
and elaborate the hardware implementation of the scheme.

###  III-A Foundation

We will consider the case that we want to perform operations (addition or
multiplication) within a ring of positive integers modulo $q$, which from here
on we denote as $\mathbb{Z}_{q}$. Let $a,b\in\mathbb{Z}_{q}$. We define
$c\in\mathbb{Z}_{q}$ as the result of ($a+b$) or ($a\cdot b$) after reduction
(mod $q$) and $\tilde{c}$ as the operation result before reduction. For the
case of addition, we can see in (3), that the maximum number we will obtain
from adding $a$ and $b$ is smaller than $2q$. This can be seen by adding the
two largest possible elements in $\mathbb{Z}_{q}$ being $q-1$. We can
therefore simply compute $(c\leftarrow\tilde{c}\mod q)$ from $\tilde{c}$ by a
conditional subtraction of $q$, as in (4).

$\forall a,b\in\mathbb{Z}_{q}:\max{\tilde{c}}=\max(a+b)=2(q-1)$ (3)
$c=(\tilde{c}\mod{q})=\begin{cases}\tilde{c}&\tilde{c}<q\\\
\tilde{c}-q&\tilde{c}\geq q\end{cases}$ (4)

For multiplication (5), the reduction becomes significantly more complex to
compute as our range of results becomes significantly larger. If we wanted to
reduce $\tilde{c}$ to $c$ we would have $q-1$ cases to consider, rendering a
solution like in (4) impractical.

$\forall a,b\in\mathbb{Z}_{q}:\max{\tilde{c}}=\max{(a\cdot b)}=(q-1)^{2}$ (5)

To solve the issue of large intermediate results, we want to consider the
following property: Let $x,y\in\mathbb{Z}$. We can see from (6), that when
computing a sum (mod $q$), we can also first compute the reduction of the
summands individually, before adding the individual elements together and
reducing the resulting sum again.

$x+y\equiv(x\mod{q})+(y\mod{q})\pmod{q}$ (6)

This property can be useful for bringing the sum closer to elements in
$\mathbb{Z}_{q}$, such that the final reduction can be performed as in (4). In
such a case, as the final reduction of the sum, we only need to conditionally
subtract $q$ if the sum is $\geq q$, since the sum of two elements in
$\mathbb{Z}_{q}$ will be at most $2(q-1)$. For this property to be useful, we
additionally need the steps $(x\mod{q})$ and $(y\mod{q})$ to be easy to
perform. We will see in the following section how we can express a
multiplication result as in (5) in such a form.

###  III-B Reduction Scheme

We define $\tilde{c}$ as the intermediate result, before reduction
$(c=\tilde{c}\mod q)$, of multiplying $(a\cdot b)$ with
$a,b\in\mathbb{Z}_{q}$. We can interpret $\tilde{c}$ as a $2n\text{-bit}$
vector, where $n=\lceil\log_{2}{q}\rceil$ is the number of bits needed to
store an element in $\mathbb{Z}_{q}$. We make the difference between
interpreting $\tilde{c}$ as an integer or bit vector explicit by using bold
$\mathbf{\tilde{c}}$ for the bit vector interpretation and normal font
$\tilde{c}$, as the integer interpretation. The relationship between
$\tilde{c}$ and $\mathbf{\tilde{c}}$ is defined as follows:

$\mathbf{\tilde{c}}=[\tilde{c}_{0},\tilde{c}_{1},...,\tilde{c}_{2n-1}],\,\tilde{c}_{i}\in\\{0,1\\}$
(7) $\tilde{c}=\sum_{i=0}^{2n-1}(\tilde{c}_{i}\cdot 2^{i})$

As introduced in (6), in a sum, we can reduce summands individually, before
reducing the resulting sum. By writing $\tilde{c}$ as a bit vector
$\mathbf{\tilde{c}}$, we can express the potentially large (relative to $q$)
number $\tilde{c}$, as a sum of powers of two. When applying the idea of
reducing summands individually from (6), to $\mathbf{\tilde{c}}$, we end up
with (8):

$c=\sum_{i=0}^{2n-1}(\tilde{c}_{i}\cdot 2^{i}\mod{q})\pmod{q}$ (8)

As previously stated, for the summation property (6) to be useful, we need to
be able to easily compute $(x\mod q)$ for a summand $x$. We can observe that
for any given bit index $i$, the value of $(2^{i}\mod q)$ is constant, since
both $q$ and $2^{i}$ are constant. These values are independent of the bit-
values in $\mathbf{\tilde{c}}$, we only changed the value of how we interpret
a bit $\tilde{c}_{i}$ in $\mathbf{\tilde{c}}$ being one from $2^{i}$ to
$(2^{i}\mod q)$ when forming the integer representation from the bit vector
representation. Summing the new bit values, results in the sum of them being
congruent to the original value of $\tilde{c}$. The difference is, that the
maximum value is now much closer to $q$, since we are summing bits whose
values are elements in $\mathbb{Z}_{q}$. We can easily precompute and store
the values of $(2^{i}\mod q)$, making the condition we’ve set, that it needs
to be easy to obtain $(x\mod q)$ in (6) true. We can rewrite an individual
summand in (8) to more clearly separate the variable bit-value from the
constant part of the expression:

$(\tilde{c}_{i}\cdot
2^{i}\mod{q})=\begin{cases}(2^{i}\mod{q})&\tilde{c}_{i}=1\\\
0&\tilde{c}_{i}=0\end{cases}$ (9)

Storing the values of $(2^{i}\mod{q})$ for all $i$ in a pre-computed table
would require us to store a total $2n$ $n\text{-bit}$ values. We can further
simplify this to just storing $n+1$ values, since the $n$ least significant
bits are already elements in $\mathbb{Z}_{q}$:

$\forall i<n:(2^{i}\mod{q})=2^{i}$ (10)

We can further group the $n-1$ least significant bits to form a single summand
since their sum will always be an element in $\mathbb{Z}_{q}$:

$q<\sum_{i=0}^{n-2}2^{i}$ (11)

Combining the things discussed so far in this section, we want to compute $c$
from the bit vector $\mathbf{\tilde{c}}$ as follows in (12). The left sum is
marked gray to point out the difference, that the left sum just requires
interpreting a bit-vector in the classical way as an integer, without needing
any actual computation. On the contrary, the right sum actually involves
summing terms of $n\text{-bit}$ numbers. Note that we deliberately include
also the $(n-1)\text{th}$ bit in the left sum, allowing it to get slightly
larger than $q$ (but smaller than $2q$) which we need to (and easily can)
account for in the final (mod $q$) step.

$c=\left({\color[rgb]{.5,.5,.5}\sum_{i=0}^{n-1}(\tilde{c}_{i}\cdot
2^{i})}+\sum_{i=n}^{2n-1}(\tilde{c}_{i}\cdot(2^{i}\mod{q}))\right)\mod q$ (12)

The final thing left up for discussion is how to compute the final reduction
(mod $q$). Thanks to the term before the final reduction just being a sum of
elements in $\mathbb{Z}_{q}$, we can again look at the largest number we
obtain from summing two elements in $\mathbb{Z}_{q}$ (3). Since the largest
value that an individual summand can have is $(q-1)$ and we are summing $n$
such terms, the maximum number that can occur will therefore be $n(q-1)$.
Since we previously allowed the term of least significant bits to become
slightly larger than, $q$ we need to account for this here by letting the
maximum number we want to reduce be $(n+1)\cdot(q-1)$. Since there are only
$n+1$ cases to consider, we can relatively cheaply build a compare logic to
subtract the correct multiple of $q$ from the intermediate result to compute
the final reduction step.

Algorithm 3 LUT based reduction scheme

1:$f[i]=2^{i}\mod q$, $n=\lceil\log_{2}{q}\rceil$, $N=2n$,
$\tilde{c}\in[0,2^{N})$

2:$c=\tilde{c}\mod q$

3:$c\leftarrow\tilde{c}\mod 2^{n}$ $\triangleright$ Take lower $n$ bits

4:$c\leftarrow c+\sum_{i=n}^{N-1}(\tilde{c_{i}}\cdot f[i])$

5:while $c\geq q$ do $\triangleright$ Maximum of $n+1$ iterations

6: $c\leftarrow c-q$

7:end while

8:return c

Alg. 3 combines all the steps described in this section. It should be pointed
out that the subtraction of $q$ in a while loop is not intended to be
implemented this way, since it would neither be very efficient nor possess a
constant execution time. In the next section, we are going to discuss how we
can put this scheme efficiently into hardware and further optimize the
summation step and the final subtraction of multiples of $q$.

###  III-C Optimization for Hardware

From the method discussed in section III-B we can build a hardware structure
as shown in Fig. 1 that resembles this functionality. We introduce $N$ as the
number of bits that need to be fed through a look-up table, which is every bit
above the $n\text{th}$ bit in the input bit-vector. In the case of reduction
after the multiplication of two elements in $\mathbb{Z}_{q}$, the resulting
bit length will be $2n$, $N=2n-n=n$. The circuit in Fig. 1 will require only a
small number of pre-computed values for the first Look-up-Table (LUT) stage
($N$ $n\text{-bit}$ values). However, the implementation will require $N+1$
$n\text{-bit}$ additions to sum up the table outputs. Additionally, in the
final reduction step, which consists of subtracting the correct multiple of
$q$ to fulfil $c<q$, all multiples of $q$ from $0\cdot q\text{ to }N\cdot q$
need to be considered.

Figure 1: The basic implementation of the proposed method before hardware
optimization

For a more efficient hardware implementation, we have two main goals with the
previously introduced basic method from Fig. 1. The first one is to reduce the
number of cases needed to be considered for conditional subtraction (final
reduction step). The second one is to reduce the number of additions needed to
sum the table outputs. What we therefore do is group our previous
$1\text{-bit}$ input LUTs from Fig. 1 into larger LUTs of $k\text{-bit
inputs}$ as illustrated in Fig. 2. $k$ is a freely choosable configuration
parameter, that can be adapted to best fit a target device’s architecture
(e.g. available LUT size of an FPGA). Note that table sizes don’t need to be
symmetric, however for better demonstration of the effect of grouping tables
together, we here chose them to be all $k\text{-bit}$ input tables. By
grouping single input bit tables into larger $k\text{-bit input}$ tables, the
number of precomputed values we need to store increases by a factor $2^{k-1}$.
The advantage on the other side is that the number of additions we needed in
the second stage decreases by (almost) a factor of $k$, since fewer total LUT
outputs need to be summed. This in turn reduces the number of cases (different
multiples of $q$) that need to be considered for the final reduction step by
(almost) a factor of $k$. The trick to reducing the number of possible
subtraction cases in the final step is, that we do not simply group and sum
all table configurations. We instead also compute a modular reduction on all
table entries, turning sums of elements in $\mathbb{Z}_{q}$, which could
normally get larger than $q$, back into elements in $\mathbb{Z}_{q}$. Bringing
the new maximum of the intermediate sum closer to $q$ (from $(N+1)\cdot(q-1)$
to $(\frac{N}{k}+1)\cdot(q-1)$).

Figure 2: The optimized implementation, forming a more balanced tradeoff of
storage and arithmetic resources. Dashed lines indicate the logic encoded
inside the LUT. The values are pre-computed and stored in the LUT for all
possible $2^{k}$ input states of the table.

The following equation adapts the original equation describing the basic
variant of the scheme (12) to the hardware optimized variant shown as a
circuit in Fig. 2:

$\hat{c}=\sum_{i=0}^{n-1}(\tilde{c}_{i}\cdot
2^{i})+\sum_{i=0}^{\frac{N}{k}-1}\left(\sum_{j=0}^{k-1}(\tilde{c}_{(n+j+ki)}\cdot
2^{n+j+ki})\mod q\right)$ (13)
$c=\hat{c}-iq\,|\,iq\leq\hat{c}<(i+1)q\,\land\,i\in[0,N/k]$

###  III-D Optimizing the Final Reduction Step

To optimize the final reduction step, we need to minimize the maximum value of
the intermediate sum $\hat{c}$. The maximum value is determined by the maximum
value of the term not feed through the LUTs $(2^{n}-1)$ (the value of all bits
in this term being one) and the sum of the largest entries in the tables
$f_{i}[x]$:

$\max(\hat{c})=(2^{n}-1)+\sum_{i=0}^{\frac{N}{k}-1}\max(f_{i}[x])\approx(2^{n}-1)+q\frac{N}{k}$
(14)

In most cases, at least one combination of the input bit vector of a table
leads to a table entry that is close to the value of $q$. We can thereby
approximate the maximum of a table as $\max(f_{i}[x])\approx q$.

One thing we can do to reduce the maximum value of $\max(\hat{c})$ that is
often worth doing, is including the $n$th bit in the tables. It could
theoretically also be forwarded to $\hat{c}$ directly, as the mod operation
does not change its value ($2^{n-1}=(2^{n-1}\mod q)$). It can however make
sense to include this bit in one of the tables. Removing the most significant
bit from the term not feed through the LUTs assures that its value remains
much smaller than $q$. This often leads to having to consider one less case in
the final reduction step.

Another optimization regards how we group the input bits into tables. The case
where this is relevant is, when the maximum value of $\hat{c}$ is just above
an integer threshold (e.g. $\max(\hat{c})=3.05$, anything $<3$ would allow us
to only need to check for multiples of $q$ up to 2 instead of 3). In this
case, we would like to lower $\max(\hat{c})$ slightly below the integer
threshold to remove one case from the final reduction step. To achieve this,
we can rearrange how we group the input bit vector $\mathbf{\tilde{c}}$ into
tables. We have so far considered the bits to be grouped in their natural
order (e.g. grouping $[c_{1},c_{2},c_{3},c_{4}]$ into 2 tables of two inputs
each would be $[[c_{1},c_{2}],[c_{3},c_{4}]]$). We can often slightly affect
the maximum value $\max(\hat{c})$ by choosing a different configuration (e.g.
$[[c_{1},c_{3}],[c_{2},c_{4}]]$). We utilize here that $\max(f_{i}[x])$ is
only approximately equal to $q$. We can vary the exact value $\max(f_{i}[x])$
by varying which input bits we group together, such that the sum of
$\max(f_{i}[x])$ becomes as small as possible.

###  III-E Numeric Example

We are now going to apply the reduction scheme introduced in III-B and its
hardware optimized variant from III-C, in a simple numeric example with only
small bit vector sizes. We let the modulus $q=13$ and the corresponding number
of bits for an element in $\mathbb{Z}_{q}$ $n=4$. We want to reduce an
arbitrarily selected $2n=8$ bit long number $\tilde{c}=210$ (binary:
$11010010$) to $c=\tilde{c}\mod q=210\mod 13=2$. Since our modulus $q$
requires four bits, we leave the four lower bits ($\tilde{c}_{0}$ to
$\tilde{c}_{3}$) untouched and precompute four single input tables (table
values $0$ and $(2^{i}\mod 13)$) for the four upper bits ($\tilde{c}_{4}$ to
$\tilde{c}_{7}$):

$i$ | $2^{i}\mod 13$
---|---
$4$ | $3$
$5$ | $6$
$6$ | $12$
$7$ | $11$
TABLE I: Table for the basic reduction scheme with single input tables

To design the final reduction step, we consider what the largest number
$\hat{c}$ would be. This would, for the basic reduction scheme, be the case if
all bits in $\mathbf{\tilde{c}}$ are $1$ and their values would be summed up:

$\max{(\hat{c})}=1+2+4+8+3+6+12+11=47$ (15)
$i_{\max}=\frac{\max{\hat{c}}}{q}=\frac{47}{13}\approx 3.6$ (16)
$c=\hat{c}-iq\,|\,i\in[0,\lfloor 3.6\rfloor]$ (17)

In this case, our final reduction step would need to consider all multiples of
$q$ from $0$ to $3$. We now have determined all needed parameters for an
instance of the reduction hardware and can start reducing numbers at runtime.
We now apply $\tilde{c}=210$ to the input of the circuit. The lower four bits
will be forwarded into the intermediate sum unchanged, which will form the
first summation term $0010\rightarrow 2$. In the upper four bits, the tables
will contribute a $0$ to the sum if the input was also $0$, or its precomputed
value if the input was $1$. Since $\tilde{c}_{4}$, $\tilde{c}_{6}$ and
$\tilde{c}_{7}$ where $1$ in the upper four bits, we add their corresponding
value from Table I, which would be: $\tilde{c}_{4}\rightarrow 3$,
$\tilde{c}_{6}\rightarrow 12$ and $\tilde{c}_{7}\rightarrow 11$. The rest of
the summation terms will be zero. We therefore obtain the intermediate sum:

$\hat{c}=2+3+12+11=28$ (18)

We now check which multiple of $q$ we need to subtract in the final step. We
can then perform the final subtraction and obtain the result:

$2q\leq\hat{c}<3q\rightarrow i=2$ (19) $c=\hat{c}-(i\ \cdot q)=28-(2\ \cdot
13)=2$ (20)

We can now observe that our end result $2$ did indeed match what we expected
from the reduction.

For the hardware optimized scheme, we let $k=2$, grouping the four upper bits
and their corresponding value from Table I into two tables of two inputs each.
The Table entries would then be:

$[\tilde{c}_{i},\tilde{c}_{i+1}]$ | $f_{1}[\tilde{c}_{4},\tilde{c}_{5}]$ | $f_{2}[\tilde{c}_{6},\tilde{c}_{7}]$
---|---|---
$[0,0]$ | $0$ | $0$
$[1,0]$ | $3$ | $12$
$[0,1]$ | $6$ | $11$
$[1,1]$ | $9$ | $(12+11)\mod 13=10$
TABLE II: Example of the table entries for the optimized scheme

We can now again determine the $i_{\max}$ parameter, which stands for the
maximum multiple of $q$ that can occur in the final subtraction step. In the
optimized scheme, this is determined by the sum of the values of the lower
bits in $\tilde{c}$ that are fed into the sum directly, plus the largest entry
in every table:

$\max{(\hat{c})}=\sum_{i=0}^{3}2^{i}+\max{(f[\tilde{c}_{4},\tilde{c}_{5}])}+\max{(f[\tilde{c}_{6},\tilde{c}_{7}])}$
(21) $=1+2+4+8+9+12=36$ $i_{\max}=\frac{\max{\hat{c}}}{q}=\frac{36}{13}\approx
2.8$ (22)

Applying $\tilde{c}=210$ to the input of this circuit will result in:

$f_{1}[\tilde{c}_{4},\tilde{c}_{5}]=f_{1}[1,0]\rightarrow 3$ (23)
$f_{2}[\tilde{c}_{6},\tilde{c}_{7}]=f_{2}[1,1]\rightarrow 10$
$[\tilde{c}_{0},\tilde{c}_{1},\tilde{c}_{2},\tilde{c}_{3}]=[0,1,0,0]\rightarrow
2$ $\hat{c}=2+3+10=15$ (24) $q\leq\hat{c}<2q\rightarrow i=1$ (25)
$c=\hat{c}-(i\ \cdot q)=15-(1\ \cdot 13)=2$ (26)

We see that we also obtain the correct result $2$ from applying this method.
We can also observe that trough grouping bits together into larger LUTs, we
have reduced the possible number of cases in final reduction from 3 cases in
(16) to just 2 in (22).

##  IV Implementation Results

In this section, we implement the scheme presented previously in III-C on a
Xilinx Zynq-7000 FPGA and present synthesis results for a potential ASIC
implementation. As a modulus, we choose the moduli of the second round post
quantum cryptography key encapsulation algorithm Crystals-Kyber [3] according
to the NIST round 3 specification [12] and the modulus of post quantum
cryptography signature algorithm Crystals-Dilithium [4]. According to their
specification, the moduli for those schemes are: $q_{\mathrm{Kyber}}=3329$ and
$q_{\mathrm{Dilithium}}=8380417$. Since Kyber uses a relatively small number
as a modulus and Dilithium a rather large number, testing the two can provide
a good estimate on how the scheme will behave for different moduli bit-widths.
The elements in $\mathbb{Z}_{q}$ will take $n_{\mathrm{Kyber}}=12$ and
$n_{\mathrm{Dilithium}}=23$ bits. We here consider a case of reduction after
the multiplication of elements in $\mathbb{Z}_{q}$, therefore we reduce bit
vectors of size $2n$. We group the $N=n+1$ MSBs into two LUTs for the case of
Kyber (7 and 6 input bit configuration) and four 6-bit LUTs in the case of
Dilithium. In FPGA implementations, we choose these groupings to make good use
of the given LUT structure within a single combinatoric logic block (CLB) of
the Zynq-7000 FPGA. This is at smallest a 5-bit LUT and at most a 8-bit LUT,
which can be formed by grouping together 6 input LUTs. Since arithmetic
operations will also be performed with LUTs in FPGAs, it is more beneficial to
choose larger (within the range of the CLBs LUTs) LUT sizes rather than
performing additional arithmetic operations. In an ASIC one may prefer to
choose slightly smaller LUTs, but more arithmetic since table sizes can be
chosen arbitrarily in ASICs and arithmetic is not built from LUTs. However, we
will stick to the 6 and 7 input LUTs for the rest of this article for both
ASIC and FPGA results.

As a reference, we also implement Barrett reduction from II, to which we will
compare our results. In total, we compare 3 implementations. ’Barrett general’
corresponds to the method shown in II-A. ’Barrett optimized’ corresponds to
II-B. Finally, LUT based corresponds to the optimized scheme shown in III-C
and Fig. 2.

###  IV-A FPGA Results

Table III shows that the LUT-based method uses a slightly larger amount of
LUTs than the general Barrett reduction, but due to the scheme not requiring
any multiplications does not instantiate any DSPs. One can also see that in
the optimized Barrett reduction the LUT count goes up as the DSP count goes
down, since we changed multiplications to shift and adds, which is performed
in the LUTs. To provide a better chip area comparison on FPGAs we also
performed an implementation run with the optimized Barrett method where we do
not allow the synthesis tool to use any DSPs for the multiplications. We see
that here the LUT count even in the optimized scheme would be more than double
compared to that of the LUT based method. This should not come as a surprise
since, also multipliers/DSPs are much fewer in numbers on FPGAs compared to
LUTs, as they consume significantly more chip area per unit, which was one of
the primary motivators for the LUT based scheme in the first place. For the
larger bit width in Dilithium, the difference between schemes (optimized and
LUT based) appear to be less significant. However, this would not generalize
to arbitrary moduli of this bit width, as the modulus of Dilithium benefits
heavily from having both a multiplication constant
$m=8396807=2^{23}+2^{13}+2^{3}-1$ and a modulus $q=8380417=2^{23}-2^{13}+1$
that can be expressed by the sum of only very few shifted input terms to
perform multiplication with them. For arbitrary moduli of this bit-width, it
is therefore better to compare it to the generalized Barrett.

| Kyber | Dilithium
---|---|---
| LUTs | DSPs | LUTs | DSPs
Barrett general | 54 | 2 | 130 | 5
Barrett optimized | 74 | 1 | 255 | 0
Barrett optimized (no DSP) | 212 | 0 | - | -
LUT based method | 82 | 0 | 200 | 0
TABLE III: FPGA implementation results in terms of resource usage

###  IV-B ASIC Synthesis Results

We perform ASIC synthesis for four reduction variants using Cadence Genus. The
purpose of this is to obtain an area comparison result, which allows for
comparison in terms of chip area for implementations that both use DSPs and
LUTs in an FPGA, such as the Barrett reduction. We have not further optimized
the LUT sizes for ASIC implementation for this test and stayed with the FPGA
optimized 6 and 7 input LUT design, which may not be the best variant for
ASICs. Since the absolute chip area will depend upon the used ASIC technology
standard, we only compare the area consumption relative to each other.

| Area normalized to Barrett optimized
---|---
| Kyber | Dilithium
Barrett optimized (synthesis tool) | 1 | 1.2 (1)
Barrett optimized (manual) | 0.97 | 1.2 (1)
Barrett general (full multiplier) | 2.7 | 7.4 (6.0)
LUT based method | 0.59 | 0.78 (0.63)
TABLE IV: Area results of ASIC synthesis normalized to the area of the
automatically optimized version of the Barrett implementation for Kyber. Terms
in () are normalized within their column

Table IV shows the result of ASIC synthesis in terms of the relative required
area on a chip. The ASIC synthesis tool can perform similar optimizations to
what we did manually in the optimized scheme II-B, since synthesis is
performed for a specific constant modulus. The synthesis tool does not need to
generate full general multipliers. We therefore also explicitly performed
synthesis for the general Barrett reduction with full multipliers, this
provides an estimate on how large the Barrett reduction module can become for
arbitrary moduli of this bit width. It can be seen that the proposed LUT based
method consumes only about 10% to 20% of the chip area compared to the general
Barrett reduction with full multipliers and about 60% of the optimized Barrett
reductions in the tested case.

##  V Conclusion

In this article, we have shown a potential alternative to the commonly used
Barrett reduction for computing modular reduction in hardware. The method
shows promising results in terms of area consumption, which was reduced
compared to optimized Barrett reduction by 39% for the tested modulus and even
90% for not specifically selected arbitrary moduli. Since especially the area
consumption of the optimized Barrett reduction varies greatly for different
moduli, the ratio by which the scheme proposed here is better (or worse) than
the Barrett method, may vary significantly depending on the modulus. For the
proposed scheme itself, the area consumption will stay relatively consistent
for different moduli of the same bit length. This can already be an advantage,
since it removes the concern of good reducibility in the process of choosing a
modulus for cryptographic schemes, which often also need to be chosen for
other properties, such as being a prime number. Since all steps of the
proposed method scale linearly to the bit width of the number being reduced,
it may also be a viable scheme for reducing large numbers. Another major
advantage of the scheme is that all operations can be performed as (close to)
$n$-bit wide logic, as opposed to the operation width being determined by the
input bit width. This shortens carry signal paths, allows for faster clocking
of the design and fast reduction even for high bit widths. Most operations
(LUTs, summation of LUT outputs) can be parallelized and are well suited for
pipeling which is another useful property for hardware implementations. As a
next step, one may want to further investigate and quantify the effects of
different parameters such as the table size $k$, modulus $q$ or size of the
number being reduced on the relative performance compared to the Barrett
method as well as the total area consumption.

## References

  * [1] Manoj Kumar and Pratap Pattnaik  “Post Quantum Cryptography(PQC) - An overview: (Invited Paper)”  In _2020 IEEE High Performance Extreme Computing Conference (HPEC)_ , 2020, pp. 1–9  DOI: 10.1109/HPEC43674.2020.9286147
  * [2] Ritik Bavdekar et al.  “Post Quantum Cryptography: A Review of Techniques, Challenges and Standardizations”  In _2023 International Conference on Information Networking (ICOIN)_ , 2023, pp. 146–151  DOI: 10.1109/ICOIN56518.2023.10048976
  * [3] Joppe Bos et al.  “CRYSTALS - Kyber: A CCA-Secure Module-Lattice-Based KEM”  In _2018 IEEE European Symposium on Security and Privacy (EuroS &P)_, 2018, pp. 353–367  DOI: 10.1109/EuroSP.2018.00032
  * [4] Léo Ducas et al.  “Crystals-dilithium: A lattice-based digital signature scheme”  In _IACR Transactions on Cryptographic Hardware and Embedded Systems_ , 2018, pp. 238–268
  * [5] Paul Barrett  “Implementing the Rivest Shamir and Adleman Public Key Encryption Algorithm on a Standard Digital Signal Processor”  In _Advances in Cryptology — CRYPTO’ 86_ , Lecture Notes in Computer Science  Berlin, Heidelberg: Springer, 1987, pp. 311–323  DOI: 10.1007/3-540-47721-7˙24
  * [6] Peter L. Montgomery  “Modular multiplication without trial division”  In _Mathematics of Computation_ 44.170, 1985, pp. 519–521  DOI: 10.1090/S0025-5718-1985-0777282-X
  * [7] Zhengjun Cao, Ruizhong Wei and Xiaodong Lin  “A Fast Modular Reduction Method” https://eprint.iacr.org/2014/040, Cryptology ePrint Archive, Paper 2014/040, 2014  URL: https://eprint.iacr.org/2014/040
  * [8] Chae Hoon Lim, Hyo Sun Hwang, Yang Jae-Dong and Seo Cho-Gu  “Fast Modular Reduction With Precomputation”  URL: https://citeseerx.ist.psu.edu/document?repid=rep1&type=pdf&doi=a683d452712796ca9ab4c99dffdfd48e1e5d195d
  * [9] Mark A. Will and Ryan K.. Ko  “Computing Mod Without Mod” https://eprint.iacr.org/2014/755, Cryptology ePrint Archive, Paper 2014/755, 2014  URL: https://eprint.iacr.org/2014/755
  * [10] Amos R.  “Modular Reduction”  In _Cryptography Arithmetic: Algorithms and Hardware Architectures_ Cham: Springer International Publishing, 2020, pp. 105–141  DOI: 10.1007/978-3-030-34142-8˙4
  * [11] Utsav Banerjee, Tenzin S Ukyab and Anantha P Chandrakasan  “Sapphire: A configurable crypto-processor for post-quantum lattice-based protocols”  In _arXiv preprint arXiv:1910.07557_ , 2019
  * [12] Information Technology Laboratory Computer Security Division  “Round 3 Submissions - Post-Quantum Cryptography — CSRC — CSRC”, 2023  URL: https://csrc.nist.gov/Projects/post-quantum-cryptography/post-quantum-cryptography-standardization/round-3-submissions