# A novel in vitro device to deliver induced electromagnetic fields to cell
and tissue cultures

Rea Ravin  Celoptics, Inc., Rockville, MD, USA  Teddy X. Cai  Randall H.
Pursley  The Signal Processing and Instrumentation Section, Center for
Information Technology, National Institutes of Health, Bethesda, MD, USA
Marcial Garmendia-Cedillos  The Signal Processing and Instrumentation Section,
Center for Information Technology, National Institutes of Health, Bethesda,
MD, USA  Tom Pohida  The Signal Processing and Instrumentation Section, Center
for Information Technology, National Institutes of Health, Bethesda, MD, USA
Raisa Z. Freidlin  The Signal Processing and Instrumentation Section, Center
for Information Technology, National Institutes of Health, Bethesda, MD, USA
Herui Wang  Neuro-Oncology Branch, Center for Cancer Research, National Cancer
Institute, National Institutes of Health, Bethesda, MD, USA  Zhengping Zhuang
Neuro-Oncology Branch, Center for Cancer Research, National Cancer Institute,
National Institutes of Health, Bethesda, MD, USA  Amber J. Giles  Neuro-
Oncology Branch, Center for Cancer Research, National Cancer Institute,
National Institutes of Health, Bethesda, MD, USA  Nathan H. Williamson
Section on Quantitative Imaging and Tissue Sciences Eunice Kennedy Shriver
National Institutes of Child Health and Human Development, National Institutes
of Health, Bethesda, MD, USA  National Institute of General Medical Sciences,
National Institutes of Health, Bethesda, MD, USA  Mark R. Gilbert  Neuro-
Oncology Branch, Center for Cancer Research, National Cancer Institute,
National Institutes of Health, Bethesda, MD, USA  Peter J. Basser

###### Abstract

We have developed a novel in vitro instrument that can deliver intermediate
frequency (100 – 400 kHz), moderate intensity (up to and exceeding 6.5 V/cm
pk-pk) electric fields (EFs) to cell and tissue cultures generated using
induced electromagnetic fields (EMFs) in an air-core solenoid coil. A major
application of these EFs is as an emerging cancer treatment modality. In vitro
studies by Novocure Ltd. reported that intermediate frequency (100 – 300 kHz),
low amplitude (1 – 3 V/cm) EFs, which they called “Tumor Treating Fields
(TTFields)”, had an anti-mitotic effect on glioblastoma multiforme (GBM)
cells. The effect was found to increase with increasing EF amplitude. Despite
continued theoretical, preclinical, and clinical study, the mechanism of
action remains incompletely understood. All previous in vitro studies of
“TTFields” have used attached, capacitively coupled electrodes to deliver
alternating EFs to cell and tissue cultures. This contacting delivery method
suffers from a poorly characterized EF profile and conductive heating that
limits the duration and amplitude of the applied EFs. In contrast, our device
delivers EFs with a well-characterized radial profile in a non-contacting
manner, eliminating conductive heating and enabling thermally regulated EF
delivery. To test and demonstrate our system, we generated continuous, 200 kHz
EMF with an EF amplitude profile spanning 0 – 6.5 V/cm pk-pk and applied them
to exemplar human thyroid cell cultures for 72 hours. We observed moderate
reduction in cell density ($<$ 10%) at low EF amplitudes ($<$ 4 V/cm) and a
greater reduction in cell density of up to 25% at higher amplitudes (4 – 6.5
V/cm). Our device can be readily extended to other EF frequency and amplitude
regimes. Future studies with this device should contribute to the ongoing
debate about the efficacy and mechanism(s) of action of “TTFields” by better
isolating the effects of EFs and providing access to previously inaccessible
EF regimes.

[pagex=2,4,6,8,10,12,14,color=gray!20,angle=45,scale=2.5,xpos=0,ypos=0]ACCEPTED
MANUSCRIPT<EMAIL_ADDRESS>Article

Previous work reported that intermediate frequency (100 – 300 kHz), low
amplitude (1 – 3 V/cm) electric fields (EFs) have an anti-mitotic effect
capable of inhibiting the growth rate of cancers cell in vitro and in the
clinic, although the mechanism of action remains unclear. Existing EF delivery
systems use capacitively coupled electrodes placed in direct contact with the
specimen holder. This direct contact generates unwanted heat and limits the
amplitude and duration of EF stimulation. We have developed a novel in vitro
system capable of continuous thermal regulation and delivery of well-
characterized (200 kHz, 0 – 6.5 V/cm pk-pk), electromagnetic fields (EMFs) to
cell cultures, paving the way for improved mechanistic biophysical studies.

## Introduction

There is an extensive body of literature describing the effects of low (DC-
kHz) and high (above 1 MHz) frequency electric fields (EFs) on living cells
(Cifra2011). Very low frequency EFs ($<$ 2 kHz) are used to stimulate
excitable cells while high frequency EFs are used to heat cells and tissues.
Intermediate frequency EFs (kHz-MHz) were previously considered to have no
biological effects. However, Kirson, et al. (Kirson2004) reported that EFs at
frequencies of $100-300$ kHz and amplitudes between $1-3$ V/cm had an anti-
mitotic effect, inhibiting the proliferation rate of human and rodent tumor
cell lines while having no effect on non-cancerous cells. The inhibitory
effect was reported to be frequency-dependent with a peak effect at 200 kHz
for glioblastoma multiforme (GBM) cell lines (Kirson2004). The inhibitory
effect was also found to increase with EF amplitude (Kirson2007). At higher EF
amplitudes ($\geq$ 2 V/cm), significant ($>$ 50%) inhibition relative to
temperature-matched controls was reported (Kirson2004, Kirson2007). Kirson et
al. labelled these EFs as “Tumor Treating Fields (TTFields)”.

GBM is the most common and lethal adult brain tumor. Prognosis remains poor.
Despite a combined treatment plan of maximal safe aggressive surgery,
chemotherapy with temozolomide, and radiation, patients’ median survival is
only 14 – 16 months from the time of diagnosis; 1 in 4 patients survive more
than 2 years and fewer than 1 in 20 survive more than 5 years (Batash2017,
Gilbert2014). One significant feature of GBM is the highly invasive character
of the cancer cells. This makes complete surgical resection impossible,
leading to a recurrence rate of almost 100% (Gilbert2014, Stupp2005). Based on
preclinical results, “TTFields” entered clinical trials (Stupp2012, Stupp2017)
and were approved by the FDA for newly diagnosed GBM cases in 2015,
representing the first new GBM treatment modality in decades. “TTFields” have
also been approved for treatment of pleural mesothelioma as of 2019
(Ceresoli2019).

Despite ongoing, theoretical, preclinical, and clinical research, the
mechanism(s) of action of “TTFields” is not well understood. Abnormal cell
mitotic morphology including blebbing followed by apoptosis (Kirson2004,
Kirson2007, Gera2015) was observed when these EFs were applied to
proliferating cells, suggesting a biophysical anti-mitotic mechanism of
action. Computational models predicted that the penetration of EFs into cells
is indeed frequency-dependent, as suggested by findings by Kirson et al.
(Kirson2004). At ultra-low frequency, EFs are largely screened and do not
penetrate the cell membrane, whereas in the low to intermediate (100 – 500
kHz) frequency range, EFs are predicted to penetrate cells and could have a
biophysical effect (Tuszynski2016, Wenger2015a, Wenger2015b, Wenger2018).

In initial preclinical studies, Kirson et al. (Kirson2004, Kirson2007)
proposed two specific mechanisms by which penetrative intermediate frequency
EFs could cause mitotic arrest and apoptosis. The first mechanism is an
interaction of these EFs with the electric dipoles of tubulin dimers which
results in changes in tubulin dimer rotational dynamics, disrupting the
formation and function of the mitotic spindle (Davies2013, Kirson2004,
Bomzon2016). The exact effects of EFs on microtubule polymerization were
predicted to depend on the length of the microtubule filament as well as the
frequency and amplitude of the applied EF. Effects on septin dimers were also
proposed (Gera2015). However, more recent computational modeling suggests that
EFs in the range of 1 – 3 V/cm are orders of magnitude too small to affect the
rotational dynamics of tubulin and septin dimers (Tuszynski2016, Wenger2015b,
Wenger2018, Berkelmann2019, Li2020).

The second mechanism of action proposed by Kirson et al. (Kirson2004,
Kirson2007) was dielectrophoretic (DEP) forces resulting in electrokinetic
effects on polarizable bio-macromolecules (i.e., those capable of producing a
counter-ion cloud). During cytokinesis, a cleavage furrow appears between the
dividing daughter cells. Computational models (Wenger2015b) predict that in
and around these furrows, non-uniform EFs (i.e., EF gradients) produce DEP
forces on charged, polarizable macromolecules and organelles (Wenger2015b,
Berkelmann2019). These DEP forces could potentially cause molecules to migrate
toward or away from the furrow, impeding cell cleavage (Tuszynski2016,
Wenger2015b, Berkelmann2019). Modeling studies incorporating Stokes drag
(Li2020), however, suggests that DEP forces are too small to cause substantial
displacement over the course of telophase. Both of these proposed mechanisms
also depend on the relative orientation of the dividing cells and the applied
EF. Parallel alignment results in a maximal effect (Wenger2015b, Wenger2018,
Kirson2007).

Recent preclinical studies support additional potential mechanisms and
downstream effects, including but not limited to mitotic checkpoint inhibition
(Kessler2018), tubulin mediated conductance effects (Santelices2017),
permeability changes (Chang2018), DNA damage (Giladi2017, Jo2018), migratory
inhibition (Silginer2017), and autophagy (Kim2019, Silginer2017). More
recently, Neuhaus et. al (Neuhaus2019) suggested voltage-gated calcium
channels and Li et al. (Li2020) suggested tumor cell membrane potential as
mediating factors for the anti-proliferative effect of intermediate frequency
EFs. While several biophysical mechanisms are plausible, none have been
conclusively demonstrated.

To advance their preclinical research on the use and mechanism(s) of action of
“TTFields”, Novocure Ltd. developed an in vitro test system (Inovitro™)
designed to deliver EFs to cell cultures using two pairs of capacitively
coupled electrodes insulated by a high dielectric constant ceramic
(Kirson2007), placed orthogonal to each other inside a cell culture dish
(Giladi2015, Kirson2007). The electrode pairs are connected to a high voltage
sinusoidal waveform generator. The stimulated electrode pair is switched
periodically to maximize the orientational efficiency of EF application
(Kirson2007). While contact between the culture dish and electrodes creates a
direct means for EFs to be delivered to cell or tissue cultures, it also
results in heat conduction from the electrodes to the dish and subsequently
from the dish to the culture media, making continuous temperature control
problematic. Joule or Ohmic heating from the EFs themselves also significantly
contributes to heating of the conductive cell media.

In order to maintain a temperature of $37{}^{\circ}$C in the culture dish, EF
application with the Inovitro™ device is performed in an 18 ∘C incubator
(Porat2017). Electric current limits and on-off operation are also imposed as
necessary to maintain temperature (Porat2017, Kirson2007). A consequence of
keeping the apparatus in this refrigerated environment is a lower surrounding
water vapor pressure, which accelerates cell media evaporation. Assuming 95%
relative humidity in the incubator, a loss of about 0.2 mL of media per day
can be expected from a 35 mm diameter dish (see Supp. Sec. 1). As a result,
media replacement is required every 24 hours according to procedures for
operating the Inovitro™device (Kirson2004, Porat2017). This rapid media loss
and replacement may have inadvertent confounding effects, such as causing
time-varying hyperosmotic bath conditions. Hyperosmotic conditions can result
in nuclear lamina buckling (Finan2010), cell cycle disruption
(Radmaneshfar2013), or protein aggregation (Fragniere2019). Another issue is
that the geometry of the EF delivery device does not lend itself to easily
predicting the electric field or current density produced within the tissue or
cell culture (Wenger2015a, Wenger2018).

To overcome some of these limitations associated with delivering EFs via
contacting electrodes, here we describe a novel device that uses
electromagnetic induction to deliver EFs/electromagnetic fields (EMFs) in
vitro, like those generated by Transcranial Magnetic Simulation (TMS) devices
(Bickford1965) but in a higher frequency range. More specifically, we propose
the use of an air-core double solenoid coil in an LC resonant circuit to
induce intermediate frequency EMF in cell or tissue cultures placed within the
coil. This non-contacting method of EMF delivery thermally isolates the
current-carrying coil from the culture dishes. Temperature uniformity and
osmolarity control of the culture medium is improved. Greater power delivery
is made possible, providing potential access to high, previously inaccessible
EF amplitudes. Furthermore, the radial, linearly varying amplitude of the
induced circumferential EF within the coil allows each culture dish to be
exposed to a known and well-characterized range of EF amplitudes or “doses”.
Induced EMFs represent a promising alternative EF delivery method for in vitro
experimentation and may better elucidate the biophysical mechanisms of action
of “TTFields” than the current state-of-the-art.

## Materials and Methods

### 200 kHz induced EMF generation

Figure 1: Experimental apparatus inside incubator. (a) Rendering of dish
sleeve containing stack of seven 35 mm cell culture dishes. Vertical positions
labelled P1 - P5. A fiber optic temperature sensor is placed into the P3 dish.
(b) Inside view of modified incubator and experimental apparatus. (1) 22 cm
$\times$ 22 cm Plexiglas window. (2) Fiber optic sensor. (3) Double-wrapped
copper coil, 3 turns, ID $=5$ cm, height $=8$ cm, attached to induction
device. (4) Infrared detector. (5) Plexiglas shelf. (c) Relative position of
control stack, off-axis from the treated sleeve to avoid the $1/r$ EMF
contribution outside the coil – see Eq. (3).

EMF generation was accomplished using an industrial 10 kW induction heater
system (DP-10-400, RDO Induction LLC, Washington NJ) based on a conventional
LC resonant circuit design (with a resonant angular frequency of 250 kHz). The
induction device was connected to a copper coil built to the following
specifications: inner diameter (ID) of 5 cm, height of 8 cm, 3 turns, double-
wrapped, with vertical orientation (Fig. 1b, 3). The induction system and coil
were jacketed and connected to a water chiller (DuraChill© DCA200,
PolyScience, Niles, IL) to remove heat and control coil temperature. In order
to place the coil inside the 95% air, 5% CO2 incubator at $37{}^{\circ}$C
(MCO-18M Multigas Incubator, SANYO Electric Co. Ltd., Osaka, Japan), a 22 cm x
22 cm hole (Fig. 1b, 1) was made in the incubator to avoid Joule heating of
the metal lining via induction currents passing through the copper pipes
connecting the coil to the induction device. The hole was covered and sealed
with Plexiglas connected to the incubator. A Plexiglas shelf was placed below
the coil (Fig. 1b, 5), replacing the metal shelf of the incubator, again to
avoid inductive heating. Plastic sleeves for cell culture dishes were made in-
house on a 3D printer (Fig. 1a). Each sleeve accommodates seven 35 mm cell
culture dishes in a vertical stack, so several culture dishes can be studied
in parallel. The middle 5 dishes are situated fully within the coil and their
positions are labelled P1 – P5. The dish in the third position from the bottom
is referred to as P3, and is shown with a fiber optic temperature sensor
placed in its center (Fig. 1a). An infrared detector records the coil
temperature (Fig. 1b, 4). EMFs were delivered continuously for 72 hours (h) in
all experiments.

The expected EF profile generated within the coil by the induction device can
be calculated from basic principles of electricity and magnetism. The magnetic
field within an air-core current-carrying coil is readily approximated by
Ampere’s law,

$B_{z}(t)=\mu_{0}I(t)\left(\frac{N}{L}\right)$ (1)

where $B_{z}(t)$ is the axial component of the applied magnetic field, B,
$\mu_{0}$ is the permeability of a vacuum, $I$ is the current amplitude, $N$
is the number of coil turns, and $L$ is the length of the coil. The direction
of the magnetic field is normal to the direction of the current and can be
determined by the right hand rule as being along the coil axis. Faraday’s law
in integral form relates the magnetic field to the electromotive force,
$\varepsilon$,

$\varepsilon=\oint_{l}\textbf{E}\cdot
d\textbf{l}=-\frac{d}{dt}\iint_{S}\textbf{B}\cdot d\textbf{S},$ (2)

where E is the induced EF, S denotes some bounded surface, and $d\textbf{l}$
is an infinitesimal arc length. If the surface is chosen to be a circular
cross section of the coil with radius $r$ from the center of the coil then the
following solutions are obtained for $|\textbf{E}(t)|$,

$|\textbf{E}(r,t)|=|E_{\theta}(r,t)|=\begin{cases}\left(\dfrac{r}{2}\right)\dfrac{dB_{z}(t)}{dt}=\mu_{0}\left(\dfrac{N}{L}\right)\left(\dfrac{r}{2}\right)\left|\dfrac{dI}{dt}\right|&r<R,\\\
\\\
\left(\dfrac{R^{2}}{2r}\right)\dfrac{dB_{z}(t)}{dt}=\mu_{0}\left(\dfrac{N}{L}\right)\left(\dfrac{R^{2}}{2r}\right)\left|\dfrac{dI}{dt}\right|&r\geq
R,\end{cases}$ (3)

where $R$ is the radius of the coil, and $E_{\theta}(r,t)$ is the azimuthal
component of the electric field. By symmetry, the other two components of the
electric field vanish. The direction of the EF is anti-parallel to the current
flowing in the coils. If the applied current is sinusoidal, e.g.,
$I(t)=I_{0}\cos(\omega t)$, then the induced EF will oscillate at the same
frequency in the circumferential direction. By Eq. (3), our induction system
delivering current at 200 kHz results in a circumferential 200 kHz EF profile
within the coil with a linearly increasing amplitude from the center of the
coil. Importantly, the radial profile of the generated EFs makes each treated
dish its own “dose titration” experiment, with little to no EFs dose or
amplitude in the center of the dish and high dose at the periphery. The
effects of EF amplitude can be isolated from any dish-to-dish experimental
confounds by analyzing radially dependent differences within a single treated
dish.

Note, the electric field outside of the coil is expected to decay by $1/r$
according to Eq. (3). The EF in the incubator outside the coil is therefore
not 0. Because of this decaying EF, the control dish sleeve is placed in the
incubator at a distance and vertical height offset (off axis from the coil)
such that little to no EF should be present, as shown in Fig. 1c.

### Temperature monitoring and control

Joule heating remains a major hurdle to delivering increased EF amplitudes
while maintaining controlled temperatures, even in the absence of conductive
heating from contacting electrodes. Heating of an electrically conductive
material by a spatially-varying EF is described by the power delivered per
unit volume

$\begin{gathered}\frac{dP}{dV}=\textbf{J}\cdot\textbf{E}=\sigma|\textbf{E}|^{2}\\\
\mathrm{by}\;\;\;\textbf{J}=\sigma\textbf{E}\;\;\;\;\mathrm{(microscopic\;Ohm^{\prime}s\;Law)},\end{gathered}$
(4)

where $P$ is power, $V$ is volume, J is current density, E is the EF, and
$\sigma$ is the conductivity of the material. In the case of a solenoid coil
(cylindrical geometry) surrounding some uniform conductive material, the rate
of heat generation, $\dot{q}$, delivered by the current-carrying coil to the
material of height, $L_{z}$, can be calculated from the within-coil part of
Eq. (3) and Eq. (4) as

$\dot{q}\equiv
P_{\mathrm{total}}=\sigma\iiint_{V}|\textbf{E}(r)|_{\mathrm{rms}}^{2}\,dV=\sigma\int_{0}^{L_{z}}dz\int_{0}^{2\pi}d\theta\int_{0}^{R}r|\textbf{E}(r)|_{\mathrm{rms}}^{2}\,dr=\frac{\pi\sigma
L_{z}m^{2}R^{4}}{16},$ (5)

where $\iiint_{V}dV$ is a volume integral and $m$ is the radial slope of the
within-coil part of Eq. (3) using peak-to-peak amplitudes such that
$|\textbf{E}(r)|_{\mathrm{rms}}=mr/(2\sqrt{2})$,

$m:=\frac{\mu_{0}}{2}\left(\frac{N}{L}\right)\left|\frac{dI}{dt}\right|_{\mathrm{pk-
pk}}\;[=]\;\mathrm{V}\cdot\mathrm{m}^{-2},$ (6)

where rms is the root-mean-square. According to Eqs. (4), (5), and (6), Joule
heat generation increases with the square of the delivered EF amplitude. Thus,
delivering higher EF amplitudes while maintaining well-controlled temperatures
requires careful regulation of environmental conditions. Not only does the
current-carrying solenoid coil require cooling, the cell cultures themselves
require negative heat flux to the surroundings to balance Joule heating. This
negative heat flux is attained by cooling the coil below the incubator
temperature. A complex interplay of incubator temperature, a cooled coil
temperature, and various thermal resistances (Fig. S2, Table S2) determines
the steady-state temperature profile within the cell culture dishes. See Supp.
Sec. 2 for a more complete discussion of heat transfer and expected
temperature profiles (Figs. S3 – S5).

A data acquisition and analysis system was used to monitor the environmental
conditions within the incubator during experiments. The system can acquire
data from as many as four temperature probes and two electric field sensors.
The data acquisition software was developed using LabVIEW (National
Instruments, Austin, TX). During operation, the system collected temperature
data from Opsens multi-channel signal conditioners (TempSens
TMS-G2-1–100ST-M1, Opsens, Quebec, Canada). Data from sensors was measured
every second and transmitted to the system computer via an RS-232 interface.
An infrared detector was used to monitor the coil temperature. Fiber optic
sensors (OTG-A-10-62ST-1, Opsens, Quebec, Canada), rather than traditional
thermistors, were used to measure temperature within the cell media to avoid
both inductive heating effects and any parasitic inductance and capacitance
potentially contributed by thermistors. Fiber optic sensors were attached to
dishes (e.g., Fig. 1b, 2) using bone wax (Lukens™ #901, Surgical Specialties,
Westwood, MA).

The water chiller set-point and temperature were also monitored and
controlled. The water temperature of the chiller was measured once per second
and transmitted to the system computer using an RS-232 interface. The
temperature set-point of the chiller was managed by the controller software
and was adjusted based on the fiber optic temperature measurements from the
center of the P3 control dish. When the measured temperature deviated from the
desired range, an alarm was triggered and the set-point temperature of the
water chiller was adjusted remotely via open-loop control. The status of the
system was monitored remotely using a cloud service (Blynk Inc., New York,
NY). The LabVIEW software communicated with the Blynk cloud service to
transmit real-time data. A Blynk cellphone app was used to receive data from
the Blynk cloud and provide the status of the ongoing experiment. The app also
had the ability to transmit data through the Blynk cloud to remotely adjust
the chiller set-point in real-time in order to maintain stable temperature
conditions throughout 72 h experiments.

### Osmolarity regulation

Cooler air near the dishes can result in accelerated loss of media due to
water vapor pressure gradients (see Supp. Sec. 1). As a result, osmolarity can
change over the course of experiments. We opted to take two preparatory steps
to mitigate osmolarity changes, as opposed to replacing cell media during the
course of the experiment. First, we used 3 mL of media rather than the
conventional 2 mL; increasing the total volume decreases the proportion of
media lost. Second, we shaped the lids of all dishes by heating the lids with
a heat gun, placing them on a sharpened ring mold, and depressing them while
malleable with a steel ball. The end result was a curved or domed lid that is
partially submerged in the media. This lid shape mitigated media loss
primarily by decreasing the exposed liquid surface area and thus the area
undergoing mass transport with the surrounding air (Fig. S1). The lid also
directs condensate back into the media volume. These changes reduced the
predicted osmolarity change by over 60% (Table S1). Note, increased media
volume and decreased exposed surface area should not affect dissolved gas (O2,
CO2, etc.) concentrations in media. According to Henry’s law, a dissolved gas
in liquid is only a function of the partial pressure of the gas above the
liquid and temperature, which are both unchanged by these modifications.
Osmolarity measurements with a Vapro 5520 vapor pressure osmometer (Wescor,
Logan, UT) were used to experimentally validate the efficacy of these
regulation measures.

### Human thyroid cell culture

We used a human thyroid cell line (Nthy-ori 3-1 Sigma, Sigma-Aldrich, St.
Louis, MO) as an exemplar test cell culture system. Cells were infected with
lentivirus green fluorescent protein (GFP) and selected for high expression of
GFP over several cycles. Cells were maintained in 75ml flasks in DMEM media
(Gibco, Thermo Fisher Scientific, Waltham, MA) supplemented with 10% FBS
(Gibco), 1% Penicillin Streptomycin (Gibco). GFP-expressing cells were
passaged twice a week and cell cultures were maintained in a 95% air, 5% CO2
incubator at 37∘C. Before each 72 h experiment, cells were detached from the
culture vessel with Accutase© (Sigma-Aldrich) and re-suspended in media at a
density of 3.33$\times$104 cells/ml. 3 ml of cell suspension was placed in
each 35 mm culture dish (Corning, Falcon) for a plating density of
1$\times$105 cells per dish. Dishes were previously coated with Poly-D-Lysine
to enhance attachment. Dishes were incubated for 24 h to allow for attachment
prior to 72 h experiments. Markings were made on the external part of the
dishes – consistent relative to manufacturer markings – to keep the
orientation of dishes consistent during treatment and imaging.

A series of cell counting calibration experiments using different initial
seeding densities ranging from 3$\times$104 to 1.2$\times$105 cells per dish
were also performed. Cells were again detached with Accutase© , plated at
different densities, incubated for 24 h, and then incubated for an additional
72 h. At the end of 72 h, the cell culture dishes were scanned using confocal
imaging, described below, and all cells were subsequently lifted and counted
using a Countess™ II FL Automated Cell Counter (Thermo Fisher). Final cell
densities ranged from $\sim 100$ to $1200$ cells/sq. mm. 3 series of
calibration experiments were performed for a total of 56 cell counts and
images. Results from these calibration experiments were used to relate imaging
observables to the experimental cell counts.

###  In situ confocal imaging and ethidium homodimer-1 and propidium iodide
staining

A unique challenge of the study design is the presence of a radial dependence
in the experimental conditions. This radial dependence precluded more obvious
approaches like lifting and counting all cells and informed our choice to use
GFP-expressing cells. Entire dishes were imaged in situ to retain radial
information, i.e., what EF amplitude treated cells experienced during the 72 h
EMF application. Stitched images were acquired and analyzed.

At the end of experiments, sleeves with dishes were removed from the coil and
dishes were gently removed from the sleeves so as not detach cells.
Orientation information was retained using dish markings. All dishes were
incubated until imaged. The dishes were mounted to the stage of a Zeiss LSM
780 laser scanning confocal microscope system. For imaging, an EC Plan-
NeoFluor 5x, 0.16 N.A. was used. All the visible parts of the 35 mm culture
dish were imaged using stitching mode, stitching 484 areas of interest with 5%
overlap. The final image size was $10,726\times 10,726$ pixels (px). For GFP,
a 488 nm excitation wavelength, a beam splitter of 488/562 nm, and emission
filters of 500 - 555 nm were used. Special care was taken to use non-
saturating parameters. Imaging parameters were kept constant across dishes and
experiments.

In 4 experiments, additional ethidium homodimer-1 (E1169, Thermo Fisher)
(EthD-1) and propidium iodide (V13241, Thermo Fisher) (PI) staining was
performed following GFP imaging. For PI, at the end of EF application and 5
minutes prior to microscope imaging, 750 $\mu$L of DMEM was removed from the
dishes and mixed with 0.5 $\mu$L of 1 mg/mL stock solution and added back to
the dish to a final concentration of $2\times 10^{-4}$ mg/mL. This procedure
was used to avoid agitating all media in the dish and losing spatial
information on floating, or incompletely adhered cells. For EthD-1, an
identical procedure was performed but for a final concentration of $0.4$
$\mu$M. No additional washes were performed to avoid loss of the floating cell
population. For both PI and EthD-1, a 561 nm excitation wavelength and
emission filters of 578 -739 nm were used. Stitching was identical.

### Image processing

Two populations were visible in images of the GFP-expressing thyroid cells: a
brighter, high fluorescence population, which might correspond to cellular
debris, mitotic cells, pre-apoptotic cells, or a rounded cell morphology; and
a dim, lower fluorescence population that corresponds to attached cells. The
primary population of interest is the attached cell population. Therefore, one
goal of image processing is to segment these bright and dim px populations for
downstream analyses. This segmentation was performed in 3 steps: contrast
enhancement and denoising (e.g., Fig. 2a and b, Fig. S7), 3-level intensity
thresholding (Figs. 2c, Fig. S8), and classification error correction (Fig.
2d, Fig. S9), outlined below. EthD-1 and PI staining images underwent a
simpler procedure of contrast enhancement and denoising (Figs. 2e and f)
followed by 2-level intensity thresholding (Fig. 2g) and size thresholding
(Fig. 2h) to identify stained dead cells. See Supp. Sec. 3 for more details
and MATLAB code. After processing, images were binned radially (Fig. 2i) to
encode radial information in the analysis. Edge regions and dish regions
occluded by fitments were excluded. The image processing pipeline is
summarized in Fig. 2.

Figure 2: Overview of image processing pipeline. (a \- d) Processing of GFP
images. (a) Raw image. (b) Contrast enhanced image. (c) Segmented image using
3-level intensity thresholding prior to corrections. (cyan arrow) Example of
“spot” classification error. (green arrow) Example of “halo” classification
error. (d) Segmented image after corrections, demonstrating good agreement
with qualitative segmentation of (a). (e \- h) Processing of PI and EthD-1
images. (e) Raw image. (f) Contrast enhanced image. (g) Segmented image using
2-level intensity thresholding. (h) Counted particles or cells based on size
thresholding. (i) Schematic representation of radial binning and exclusion of
occluded regions on an example stitched dish image. Shaded region is included.
Image has been linearly contrast enhanced ([0.01, 0.2] $\rightarrow$ [0.1, 1])
and down sampled to 20% resolution to aid visualization.

#### Contrast enhancement and denoising

The contrast enhancement step (Figs. 2b and f, Fig. S7, Supp. Sec. 3.1) was
based on a morphological transform method (Hassanpour2015, Soille2004) (Supp.
Sec. 3). An image, $A$, can be contrast enhanced by computing

$A_{\mathrm{Enhanced}}=A+\left[A-(A\circ B)\right]-\left[(A\bullet
B)-A\right]$ (7)

where $B$ is a structuring element, $\circ$ denotes an opening operation, and
$\bullet$ denotes a closing operation. Opening is erosion followed by dilation
while closing is dilation followed by erosion. The middle term on the right-
hand side of Eq. (7) is a ‘top hat’ of the original image, weighted to
brighter regions, and the rightmost term is a ‘bottom hat’, weighted to darker
regions. Morphological contrast enhancement can therefore be conceptualized as
‘adding’ to bright regions, and ‘subtracting’ from dark regions
(Hassanpour2015). All images underwent morphological contrast enhancement
using a ‘disk’ shaped structuring element with a 3 px radius. Operations were
performed using MATLAB 2019a (MathWorks, Natick, MA) and the imtophat() and
imbottomhat() functions found in the Image Processing Toolbox. After contrast
enhancement, a denoising step was performed with a 3$\times$3 Wiener filter.
Identical contrast enhancement and denoising was performed for all GFP,
EthD-1, and PI images.

#### Segmentation via thresholding

To segment the contrast enhanced GFP image, two intensity thresholds were
needed: a lower threshold to separate the background and dim px population and
an upper threshold to separate the bright and dim px populations (Supp. Sec.
3.2). The upper threshold (0.3727) was selected by applying Otsu’s method
(Otsu1979) with 3 clusters to intensity histograms drawn from a central region
of all control dish images from positions P2, P4, and P5 and then taking the
average of the larger of the two returned intensity thresholds. The lower
intensity threshold (0.0635) was selected by first filtering the contrast
enhanced and denoised image region using a $3\times 3$ matrix of ones, i.e., a
standard deviation filter. The filtered image was then thresholded using
Otsu’s method with 2 clusters. The px population with intermediate standard
deviation values was identified as corresponding to attached cell edges (Fig.
S8a). An intensity threshold was then selected such that no more than 25% of
the identified edges were excluded by the threshold (Fig. S8c). This process
was repeated for all central regions of P2, P4, and P5 control dish images.
The averages of the selected lower and upper intensity thresholds (Fig. S10)
were used for all images (e.g., Fig. 2c, Fig. S9).

To segment the contrast enhanced EthD-1 and PI images, Otsu’s method was used
to generate a single intensity threshold. After quantizing (see Fig. 2g),
contiguous regions larger than 3 px and smaller than 50 px were counted as
stained cells (Fig. 2h).

#### Classification error correction

After applying thresholds and quantizing filters to cluster the GFP image into
background, dim, and bright populations, two types of classification errors
were corrected (Supp. Sec. 3.3, Fig. S9). The first was a “halo” effect around
bright regions caused by intensity drop-off such that the edges of bright
regions were incorrectly classified as dim px. The second was a “spot” effect
in the center of dim regions caused by intensity peaks. To correct the first
type of error, all bright regions exceeding 8 px in size with a surrounding px
population of at least 60% dim px were expanded by one px layer to fill the
surrounding dim px (Figs. 2c and d, green arrow), i.e., a flood fill. To
correct the second type of error, all small bright regions ($<$ 8 px) that
were fully enclosed by dim px were converted to dim px (Figs. 2c and d, cyan
arrow). Two passes of these corrections were performed.

#### Radial binning and exclusion

Data was analyzed as a function of radius by constructing radial bins or
bands. For GFP images, an arbitrary number of equal annular width bands, 25,
was used to divide the total image area starting from the dish center. The
choice of band number results in approximately 0.72 mm wide bands. The bottom
portion of the dish is occluded by a fitment and also shows irregularities. To
avoid this area, only an upper portion of the band was kept from band 17 and
on. The edge of the dish was also excluded due to irregularities and observed
temperature differences: bands 23 – 25 were discarded from all analyses (Fig.
2i). For the EthD-1 and PI stained images, a similar binning procedure with 12
bands of approximately 1.46 mm annular wide was used, with exclusion of lower
regions after band 8, and a complete exclusion from band 10 onward.

### Study design

In each experiment, 10 dishes were simultaneously seeded. Dishes were then
randomized to either the control or treated dish sleeve/stack. Pairs of dishes
were assigned to the same height position (P1 – P5). For the dishes assigned
to P3, one dish is discarded and a dish with only media is used to
continuously probe temperature in the treated dish sleeve. After 24 h
incubation followed by 72 h of incubation and continuous 200 kHz induced EMF
application to the treated dish sleeve, dishes were imaged and processed as
described. Using data from cell counting calibration experiments, the imaged
dim px density was related to cell density using an empirical fit, shown
later. Each treated dish’s radially dependent cell density was normalized to
its position matched (e.g., P2 treated vs. P2 control) control dish. In 4
experiments, additional PI and EthD-1 staining, imaging, and image processing
was performed to clarify the role of cell death in results.

Each experiment generated 3 replicates of dish pairs at P2, P4, and P5. Pairs
of dishes assigned to P1 were discarded from all analyses due to a substantial
temperature difference, discussed later. Each dish pair was treated as a
single experimental observation of the effect of 200 kHz EFs. Data from Kirson
et al. (Kirson2007) suggest that treatment of proliferating human cells with
an EF of 200 kHz and 1 V/cm results in a decrease of roughly $20\pm 10\%$ in
cell number after 72 h compared to untreated cells. Presuming an effect of
this magnitude, 18 replicates are needed to obtain a conservative statistical
power of $0.95$ and a confidence level of $\alpha=0.01$. 12 experiments were
performed using our system, resulting in $N=36$ replicates. In 2 of these
experiments (6 replicates), PI staining and imaging was performed after
primary imaging. In 2 other experiments, EthD-1 staining and imaging was
performed for a total of $n=$ 12 out of 36 replicates assessing cell death.

### Statistical analysis and methods

All statistical analyses were performed in R 3.6.2. The linear fit for the
experimentally measured EF amplitude was performed using the typical linear
least squares (lm) approach. The bounded exponential fit for dim px density
vs. cell density calibration (of the form: $[\mathrm{Dim\;Px\;Density}]\sim
1-\exp{\\{\beta_{1}\times[\mathrm{Cell\;Density}]+\beta_{0}\\}}$) was
performed using the Gauss-Newton nonlinear least squares (nls) approach.
Initial guesses of $\beta_{1}=-0.001$ and $\beta_{2}=0$ were provided.
Bootstrapping for 95% confidence intervals was performed using the bias-
corrected and accelerated (BCa) approach provided in the R boot package (ver.
1.3-24). For bootstrapping, 1000 re-samplings were performed, following
convention. Two-sample independent $t$-tests were performed using the base R
t.test function assuming unequal variances, i.e., Welch’s $t$-test.

## Results

### Delivery of 200 kHz EFs up to $\sim$ 6.5 V/cm pk-pk to cell cultures via
induced EMFs

Figure 3: Measurements of EF amplitude as a function of pickup coil
radius/radial distance at different dish stack positions. Black line is a
zero-intercept linear least squares fit using data from all positions ([EF
Amplitude] = 4.242 $\times$ [Radial Distance], Std. Error = 0.071, Adj. R2 =
0.9902, p $<$ 0.001).

The EF amplitude and frequency delivered to the treated dish sleeve were
measured using pickup coils (i.e., loops) of several different radii composed
of insulated twisted wire pairs to minimize stray inductance. These pickup
loops were hot glued to dishes, aligned centrally, and placed within the
sleeve. The generated peak-to-peak voltage (Vpk-pk) was measured with a two-
channel oscilloscope (Hantek6022BL, Hantek Electronic Co., Ltd., Qingdao,
China) and EF amplitude was determined by $|\textbf{E}|=\mathrm{V_{pk-
pk}}/(2\pi r)$ where $r$ is the loop radius. Results from these measurements
are shown in Fig. 3, and agree with the expected linear EF amplitude profile
predicted within the coil from Eq. (3). 200 kHz oscillation was also observed
in all measurements, confirming the linearity of the system. Note, EF
amplitudes throughout refer to pk-pk values. We obtained a maximal EF
amplitude of $\sim$ 6 – 7 V/cm at the periphery of the dish with a power level
of 3.19 kW. The power level is a consequence of 220 VAC operation and the
inductive load of the culture dishes and fixture system. We report a method of
delivering EFs with a large dynamic amplitude range and well-characterized
amplitude profile.

### Temperature and osmolarity control during 72 h continuous 200 kHz EMF
delivery

Figure 4: Quantification of temperature stability and homogeneity. (a)
Temperature log data for the entire duration of a single 72 h experiment from
(top) the chiller water loop, (middle) the infrared detector reading of the
coil, and (bottom) the fiber optic sensor placed in the center of the dish at
P3. Dashed lines indicate $37\pm 0.5{}^{\circ}$C. Note the fluctuations in
coil temperature and chiller set-point adjustment throughout the experiment.
(b) Adjusted temperature data from 3 fiber optic probes placed in different
positions of the dish at P3 during the first 17.5 h of an experiment. (c)
Temperature data from fiber optic probes placed in the centers of dishes P1,
P2, P4, and P5, expressed as a difference from a simultaneous reading of a
probe in the center of the dish P3 with dashed lines indicating $0\pm
1{}^{\circ}$C.

These cell cultures are typically grown at a temperature of $37{}^{\circ}$C.
Extensive steps were taken to maintain this desired temperature in the treated
dish sleeve as growth rate is known to be function of ambient temperature. For
the present experimental setup, $m=4.242$ V/cm2 (slope of the fit Fig. 3),
such that $0.1$ W $\equiv 3.6$ kWh of Joule heating is expected per dish by
Eq. (5) (see Supp. Sec. 2.4, Eq. S23). To balance this heat generation, the
coil is jacketed and actively cooled by an external water chiller to $\approx
22\pm 1{}^{\circ}$C, cooling the surrounding air to keep the dish at P3’s
temperature at a homogeneous $37\pm 0.5{}^{\circ}$C. No inadvertent effects
like Rayleigh-Bénard convection are expected using this temperature control
method (Supp. Sec. 2.5, Eq. S24, Fig. S6). This method of temperature control
was able to maintain the desired P3 dish temperature, without interruption and
without fail, for all 72 h experiments. Exemplar temperature logs for the P3
fiber optic sensor, the infrared coil recording, and the chiller loop over the
course of one experiment are shown in Fig. 4a.

To assess the homogeneity of temperature during a stimulation experiment,
measurements from different locations within the P3 dish (Fig. 4b) (within-
dish variability) and measurements from the center of different dish positions
(Fig. 4c) (between-dish variability) were taken. For the measurements between
different dish positions, data is presented as a difference from a
simultaneous P3 dish measurement. In Fig. 4b, stimulation was started at time
$=0$. Temperature moves in tandem, demonstrating no time-lag in temperature
adjustment. Readings show that temperature within the dish is similar, but the
edge of the dish is $\approx$ 0.5∘C cooler. This was likely the result of
convective heat transfer with the cooler surrounding air (see Supp. Secs. 2.1
- 2.4). The edges of all dishes were excluded in downstream analyses as
described in Fig. 2i to ensure the integrity of within-plate information. In
Fig. 4c, data is shown from segments of experiments that are underway and at
steady state thermal conditions. Temperature measurements from the centers of
P2, P4, and P5 remained within 1∘C of measurements from the center of P3 and
were stable. Measurements from the center of P1, however, were several ∘C
cooler and so P1 was excluded from all analyses. Overall, temperature control
was satisfactory and allowed us to retain most data with confidence that
temperature differences were minimal and unlikely to contribute to observed
differences between experimental conditions.

Figure 5: Quantification of osmolarity regulation improvements. (a) Lid
shaping process with steel ball impressing a pre-heated lid placed on a
sharpened ring mold. Note exposed surface area reduction in shaped lid. (b)
Osmolarity changes after a 72 h experiment in different conditions and
configurations with 3 repetitions. Error bars = $\pm 1$ SD.

Osmolarity was measured after a 72 h experiment for control and treated dishes
in 3 different dish configurations: (1) flat lid, 2 mL media, (2) flat lid, 3
mL media, and (3) curved lid, 3 mL media, to quantify the efficacy of the
described mitigation methods (Fig. 5, Table S1). 3 repetitions of each
condition (P2, P4, P5 from one experiment) were performed. No cells were
present in the media during these osmolarity experiments. Results in Fig. 5b
are expressed as a percentage change from the measured initial average
osmolarity of 314 mOsm. The combined mitigation techniques reduced the
osmolarity increase in the treated condition from about 39% to 17% and in the
control condition from about 10% to 6.5% (Supp. Sec. 1, Eq. S10). The
reduction in osmolarity increase translated to a final osmolarity of about 370
mOsm in the treated condition compared to about 335 mOsm for the control
condition. The mitigation was imperfect but nonetheless a marked improvement
from a final osmolarity of about 440 mOsm in the unmitigated, treated
condition.

### Estimating counted cell density using segmented px density obtained via in
situ imaging of GFP-expressing thyroid cells

Figure 6: Translating dim px density to cell density. (a) Data from
calibration experiments. Dim px density plotted vs. cell density determined
from cell counts. Bounded exponential fits of individual experiments shown in
colored lines. A fit with pooled data ($n=56$) is shown in black
($[\mathrm{Dim\;Px\;Density}]=1-\exp\\{-0.00109\times[\mathrm{Cell\;Density}]-0.0784\\}$,
Residual sum-of-squares = 0.1922, Spearman’s $\rho$ = 0.949, p $<$ 0.001). (b)
Distribution of dim px densities in included radial bands of control dishes
($n=36\times 23$). The mean value for the distribution is shown with
corresponding predicted cell density using the pooled fit.

GFP images from calibration experiments were processed as described in Fig. 2.
Data from these calibration experiments is expressed as total (whole dish) dim
px density (px/px) vs. cell count per area (cells/sq. mm) in Fig. 6a. Grouped
data corresponded well to a bounded exponential fit, allowing for the robust
conversion of px densities to approximate cell densities. The bounded
exponential form of the fit may arise due to changes in the size and
morphology of thyroid cells as they approach confluence. To demonstrate that
our ‘calibration curve’ was valid for the range of data observed, a
probability density distribution of dim px densities in all included radial
bands of all P2, P4, and P5 control dishes is shown alongside the fit (Fig.
6b). Observed dim px densities fell within the range of the calibration curve.
The mean value of this distribution is also indicated. The mean dim px density
is 0.374 and the corresponding cell density from the fit is 358 cells/sq. mm.

### Application of 200 kHz induced EMF for 72 h causes moderate EF amplitude-
dependent reduction in human thyroid cell density, especially at EF amplitudes
exceeding 4 V/cm

Figure 7: Example of treated vs. control dish and data processing and
representation. (a – c) Example of (a) control dish and zoomed, processed
regions near the center (b, yellow) and near the periphery (c, magenta). Whole
dish image is adjusted as in Fig. 2i for visualization. (d – f) Corresponding
treated dish with the same zoomed regions. An overlay of the applied 200 kHz
EMF amplitude and azimuthal direction is shown. (g – h) Radially binned data
comparing control and treated dishes shown in parts (a - f). (g) Raw dim px
density converted to cell density (h, left axis) via the fit in Fig. 6, and
subsequently normalized to the mean of the control curve (h, right axis).

A pair of control and treated culture dishes is shown in Figs. 7a – b. The
displayed images provide an example of a radially dependent cell density
reduction in the treated dish. Both a central and peripheral region of the
control and treated dishes are shown at 100% resolution (Figs. 7b – c, e – f)
after the described image processing pipeline (Figs. 2a – d) to demonstrate
this radially dependent reduction in dim px density and thus cell density.
Conversion from radially binned dim px density to a normalized cell density is
also shown in Figs. 7g – h using the fit in Fig. 6.

Figure 8: Summary of treatment results presented as normalized cell density
curves converted from dim px densities. (a) Control vs. treated curve
($N=36$). Ribbons = 95% CI determined by percentile bootstrapping method with
1000 re-samplings. Point-wise t-tests: $**=$ p $<0.01$, $***=$ p $<0.001$.
Clinically relevant “TTFields” range shown for comparison. (b) Difference
between point-wise control vs. treated normalized cell densities bootstrapped
in the same way. Dashed line at 4 V/cm indicates apparent regime change.

Data from all replicates ($N=36$) was converted to normalized control vs.
treated cell density curves as shown in Fig. 7h. We address these normalized
dim px density results first, and bright px density and EthD-1 and PI staining
later. Average curves with bootstrapped 95% confidence intervals (CI) are
shown in Fig. 8a. The bootstrapped difference between curves is shown in Fig.
8b. Results demonstrate that significant differences in cell density are
present throughout the dish, even at EF amplitudes thought to be sub-
therapeutic ($<$ 1 V/cm) (Kirson2007). These differences near the center of
the dish cannot be clearly attributed to EF stimulation and could potentially
be the result of conditions present throughout the dish, such as the final
osmolarity difference between control and treated conditions discussed in Fig.
5, mild temperature gradients (see Fig. 4, Supp. Sec. 2), or the alternating
magnetic field (see Discussion). At higher EF amplitudes ($>$ 4 V/cm) an EF
amplitude-dependent effect emerges and the difference between the control and
treated conditions increases linearly with EF amplitude from about 10% to 25%
from 4 to 6.5 V/cm. At lower EF amplitudes in the same dish, a constant
difference is instead observed, suggesting the presence of heterogeneous
effects in different amplitude regimes. Note, also, that this is a magnitude
cell density decrease in the treated condition, suggesting that the peripheral
decrease is not a result of cell motility. A motility effect would instead
cause a relative increase in central density. Overall, these results suggest a
modest effect on human thyroid cell line density by 200 kHz EMF, and,
consistent with the literature, the effect is more pronounced at higher EF
amplitudes.

### Reduction in human thyroid cell density is not sufficiently explained by
cell death

Figure 9: Cell death and bright px density in response to treatment. (a \- b)
Dead cell density curves and differences constructed from pooled EthD-1
staining ($n=6$) and PI staining ($n=6$) results processed as described and
expressed as a percent of control mean cell density (358 cells/sq. mm).
Ribbons = 95% CI. No point-wise two-sample independent $t$-tests were
significant. (c \- d) Bright cell density curves and difference ($n=36$),
expressed as percent of control mean dim px density (0.374 px/px). Again, no
point-wise two-sample independent $t$-tests were significant.

The observed reduction in cell density can be explained by a reduction in
proliferation, an increase in cell death, or some mixture of both mechanisms.
To explore their relative roles, bright px density and dead cell density were
both quantified, shown in Fig. 9. Bright px density, which may correspond to
cellular debris, pre-apoptotic cells, mitotic cells, or a rounded cell
morphology, was quantified in the same way as dim px density (Fig. 2d) but
considering the brighter (violet) segmented population. Dead cell density was
quantified as described in the methods. Data from experiments with EthD-1 and
PI staining are pooled. Both bright px density (Fig. 9c – d) and dead cell
density (Fig. 9a – b) are plotted as a percentage of the control mean values
indicated in Fig. 6 to contextualize the relative contribution of bright px
and cell death to the effect observed in Fig. 8. Results indicate that while
the treated condition exhibits some qualitative differences in the radial
curves for bright px and dead cell density, these differences are not
statistically significant at the replicate number used. Moreover, as a
percentage of the control means in Fig. 6, the differences are small ($<$ 1%)
compared to the 10 – 25% differences observed in Figs. 8a and b. Bright px and
dead cells therefore cannot sufficiently explain the observed effects,
especially at higher EF amplitudes. The observed effect appears to be anti-
proliferative rather than destructive or morphological based on these
preliminary mechanistic explorations.

## Discussion

In this work, we detail the development of methods to deliver moderate
amplitude, intermediate frequency EFs using induced EMFs. We were able to
experimentally validate the delivery of a spatially varying EF profile and
quantify substantial experimental improvements in thermal and osmolar
regulation. We also report novel results of 72 h application of continuous,
200 kHz, 0 to 6.5 V/cm EMF stimulation on a non-cancerous, proliferating human
cell line. Results were analyzed using a newly developed image processing
pipeline designed specifically to characterize the in situ images necessitated
by the application of a spatially varying EF “dose” profile. Our results
suggest that (1) 200 kHz EMF have a moderate maximal effect (< 25%) on non-
cancerous cells in this amplitude range, (2) the effect of 200 kHz EMF on cell
density is likely anti-proliferative (i.e., anti-mitotic) rather than
destructive, and (3) the effect on cell density increases with increasing EF
amplitude in a higher EF regime (> 4 V/cm).

Our exploratory results corroborate and extend prior findings on the effects
of “TTFields” on proliferating, non-cancerous cell lines. Kirson et al.
(Kirson2004) reported that baby hamster kidney (BHK) fibroblasts grown in 0.1%
fetal calf serum exhibited no decrease in growth rate during 24 h application
of 1.2 V/cm pk, 100 kHz EFs. Similarly, Jo et al. (Jo2018) reported little to
no reduction in cell count for normal human skin cells (CCD-986sk) after 72 h
application of up to 1.5 V/cm pk, 200 kHz EFs. Our results agree with these
findings. Within and below the therapeutic EF amplitude range of 1 – 3 V/cm
pk-pk, we found that human thyroid cells (Nthy-ori 3-1 Sigma) did not exhibit
EF amplitude-dependent growth reduction. We then extend prior findings by
showing that, for this cell line, an EF amplitude-dependent effect is observed
at $>$ 4 V/cm pk-pk – which is well beyond the amplitudes used by Kirson et
al. (Kirson2004) and Jo et al. (Jo2018). These results suggest that there may
be a range of EF amplitudes within which some non-cancerous cells are
unaffected. Further studies may clarify whether such an EF amplitude threshold
exists for other non-cancerous cell lines such as CCD-986sk.

Our findings also provide preliminary support for anti-mitotic mechanisms of
action such as the mitotic spindle disruption reported by Giladi et al.
(Giladi2015), rather than cell death-related mechanisms (i.e., autophagy or
necroptosis). However, these findings do not necessarily generalize to
cancerous cell lines. Indeed, it is known that the effects of alternating EFs
are highly variable between cell lines (Kirson2004, Jo2018, Gera2015,
Giladi2015, Silginer2017). Future studies using cancerous cell lines are
necessary to substantiate any mechanistic claims. Nonetheless, these results
constitute a proof-of-concept for our device and lay the methodological and
procedural foundation for future studies of glioma and other cancer cell
lines. Importantly, temperature results in Fig. 4 and Fig. S3 confirm that our
device can potentially operate at even greater power (and thus EF amplitudes)
by simply lowering the coil temperature, providing the capability to
investigate the effects of previously unexplored EF amplitudes.

More specifically, our device is well-suited to begin clarifying biophysical
mechanisms that can be assessed at an “endpoint”. Some of the proposed
mechanisms can only be studied with simultaneous EF delivery and imaging
(e.g., Ref. (Giladi2015)), which is not feasible with the current
configuration of our device. Other mechanisms, however, are easily assessed at
the experimental endpoint. For example, Kim et al. (Kim2019) reported an
autophagy effect in glioblastoma cells mediated by deactivation of the
Akt2/miR29b axis. Our device can readily study radially-dependent cell death,
as reported via imaging observables here. Complementary genomic or
transcriptomic studies can also be performed by merely lifting cells from
different radial positions within the dish. Thus, we can begin to validate
findings such as those reported by Kim et al. (Kim2019) using our alternative
EF delivery platform. As a second example, reported DNA damage effects
(Giladi2017, Jo2018) can likewise be assessed by performing the appropriate
assays (e.g., an alkaline comet assay) on cells lifted from different radial
positions. We reiterate that our device has the unique advantage of each dish
representing an internally controlled EF dose titration experiment.

Although treated interchangeably with EFs throughout, it should be noted that
using induced EMFs to generate 200 kHz EFs results in an accompanying magnetic
field that was not present when using previous delivery methods. At the power
levels used, a vertically alternating uniform magnetic field of about 0.6
millitesla is expected to be produced within the coil according to Eq. (1).
This magnetic field may cause its own biological effects. Recent research
suggests that high gradient and amplitude magnetic fields can affect cell
division and resting membrane potentials (Zablotskii2016). The magnetic field
produced here, however, is orders of magnitude smaller than those shown to
have biological effects. Moreover, because biological cells and tissues have a
diamagnetic susceptibility nearly identical to that of water (Schenck1996) and
the magnetic field within a solenoid coil is known to be uniform (and not
spatially varying), we do not anticipate that the applied magnetic field will
result in any significant magnetic field gradients at interfaces. We
acknowledge, however, that magnetic field effects are one possible explanation
for the radially uniform level “offset” observed in Fig. 8 which cannot be
ruled out at this time.

Finally, while the proposed methods offer clear improvements over previous EF
delivery methods, our experimental test system is not without practical
limitations. The throughput of the system is limited by the use of a single-
coil solenoid as the EF applicator. Higher throughput experiments would be
bottle-necked by the space requirements of the current platform, which
involves a large chiller and induction heater. Another practical drawback is
the power requirement; the present 3.19 kW power requirement approaches the
battery capacities of industry-leading electric vehicles over a single day
($\approx 75$ kWh). Future developments should attempt to miniaturize the
experimental platform and introduce parallelism to enable combinatorial
experiments. Alternative coil architectures, chiller loops, multi-coil array
designs, controllable capacitor banks, etc., are all worthwhile engineering
improvements that could enhance the utility, applicability, and efficiency of
this platform. These designs could be motivated by more detailed finite
element modeling (FEM) of the various components of this test system, which is
beyond the scope of this present work.

## Conclusion

We report methods to deliver EMFs to cell and tissue cultures using
electromagnetic induction, which overcomes prior limitations associated with
the delivery of 200 kHz EFs using contacting electrodes. The proposed strategy
for inductive EMF delivery eliminates conductive heating of culture dishes and
improves temperature control and regulation. Additional design enhancements of
our culture dishes improved osmoregulation by mitigating media loss. We report
results of stimulation of human thyroid cell line cultures with 200 kHz EMFs
exceeding 3 V/cm pk-pk and as high as 6.5 V/cm. Results support a moderate
anti-proliferative effect, especially at EF amplitudes exceeding 4 V/cm. The
effect is unlikely to be due to cell death based on secondary PI and EthD-1
staining. In summary, the 200 kHz EMF delivery method reported here is a
carefully vetted and well-characterized in vitro experimental platform for
future exploration of the effect of intermediate frequency EMFs on living (and
non-living) systems.

## Author Contributions

R.R., P.J.B. conceived the study. R.R., P.J.B., T.X.C. designed the study.
R.H.P., M.G.C., T.P. contributed to electric field delivery and monitoring
setup. R.R., H.W., A.G. contributed to cell culture methods. R.R. performed
all experiments and acquired data. T.X.C. designed the analysis, analyzed
data, and wrote the supplementary material. R.Z.F. contributed to image data
analysis. T.X.C., R.R., P.J.B. wrote the manuscript. N.H.W. gave technical
advice. M.R.G. and Z.Z. provided cancer and cell biology expertise and
contributed to the discussion. P.J.B. supervised the project. All authors
edited the manuscript.

## Acknowledgments

Light microscopy was performed at the NICHD Microscopy & Imaging Core with the
assistance of Dr. Vincent Schram and Mrs. Lynne Holtzclaw. We would also like
to thank Dr. Nicole Morgan for helpful discussions. The authors declare no
conflicts of interest. Support for this study comes from the Intramural
Research Program (IRP) of the Eunice Kennedy Shriver National Institute of
Child Health and Human Development.

## Supporting Material

Supporting material accompanies this article and contains supplemental
sections concerning: (1) estimated media loss calculations, (2) heat transfer
and temperature profile calculations, (3) image processing and threshold
selection implementation. References (Welty2015, Herning1936, Schwertz1951,
Teske2005, Arnold1994, Chen2009, Chapman1990, Wilke1950, Harper2000,
Churchill1975, McAdams1985, Bird2006, CRC, Mittal2007) appear in the
supporting material. The supporting material can be found by visiting BJ
Online at http://www.biophysj.org.


Unlike ordinal-utility matching markets, which are well-developed from the viewpoint of both
theory and practice, recent insights from a computer science perspective have left
cardinal-utility matching markets in a quandary.
The celebrated pricing-based mechanism for one-sided cardinal-utility matching markets due to
Hylland and Zeckhauser [22], which had long eluded efficient algorithms, was finally
shown to be PPAD-complete [31, 13].

This led us to ask the question: is there an alternative, polynomial time, mechanism for
one-sided cardinal-utility matching markets which achieves the desirable properties of HZ, i.e.    (ex-ante) envy-freeness (EF) and Pareto-optimality (PO)?
In this paper we show:

* The problem of finding an EF+PO lottery in a one-sided cardinal-utility matching
market is PPAD-complete.

* A $(2 + \epsilon)$-approximately envy-free and (exactly) Pareto-optimal lottery can be
found in polynomial time using Nash bargaining.

We also present several results on two-sided cardinal-utility matching markets, including
non-existence of EF+PO lotteries as well as existence of justified-envy-free and weak
Pareto-optimal lotteries.

§ ONE-SIDED MATCHING MARKETS

§ TWO-SIDED MATCHING MARKETS

[1]
Atila Abdulkadiroglu and Tayfun Sönmez.
Random serial dictatorship and the core from random endowments in
house allocation problems.
Econometrica, 66(3):689–701, 1998.

[2]
Saeed Alaei, Pooya Jalaly Khalilabadi, and Eva Tardos.
Computing equilibrium in matching markets.
In Proceedings of the 2017 ACM Conference on Economics and
Computation, EC '17, page 245–261, New York, NY, USA, 2017. Association
for Computing Machinery.

[3]
Itai Ashlagi and Peng Shi.
Optimal allocation without money: An engineering approach.
Management Science, 62(4):1078–1097, 2016.

[4]
Robert J. Aumann.
Markets with a continuum of traders.
Econometrica, 32(1/2):39–50, 1964.

[5]
Haris Aziz and Ethan Brown.
Random assignment under bi-valued utilities: Analyzing
hylland-zeckhauser, nash-bargaining, and other rules.
arXiv preprint arXiv:2006.15747, 2020.

[6]
Garrett Birkhoff.
Tres observaciones sobre el algebra lineal.
Univ. Nac. Tucuman, Ser. A, 5:147–154, 1946.

[7]
Anna Bogomolnaia and Hervé Moulin.
A new solution to the random assignment problem.
J. Econ. Theory, 100:295–328, 2001.

[8]
Anna Bogomolnaia and Hervé Moulin.
Random matching under dichotomous preferences.
Econometrica, 72(1):257–279, 2004.

[9]
Eric Budish.
The combinatorial assignment problem: Approximate competitive
equilibrium from equal incomes.
Journal of Political Economy, 119(6):1061–1103, 2011.

[10]
Ioannis Caragiannis, Aris Filos-Ratsikas, Panagiotis Kanellopoulos, and Rohit
Stable fractional matchings.
In Proceedings of the 2019 ACM Conference on Economics and
Computation, EC '19, page 21–39, New York, NY, USA, 2019. Association for
Computing Machinery.

[11]
Ioannis Caragiannis, Kristoffer Arnsfelt Hansen, and Nidhi Rathi.
On the complexity of pareto-optimal and envy-free lotteries.
arXiv preprint arXiv:2307.12605, 2023.

[12]
Ioannis Caragiannis, David Kurokawa, Hervé Moulin, Ariel D. Procaccia,
Nisarg Shah, and Junxing Wang.
The unreasonable fairness of maximum nash welfare.
ACM Trans. Econ. Comput., 7(3), sep 2019.

[13]
Thomas Chen, Xi Chen, Binghui Peng, and Mihalis Yannakakis.
Computational hardness of the Hylland-Zeckhauser scheme.
In Proceedings of the 2022 Annual ACM-SIAM Symposium on Discrete
Algorithms (SODA), pages 2253–2268. SIAM, 2022.

[14]
Richard Cole and Yixin Tao.
On the existence of pareto efficient and envy-free allocations.
Journal of Economic Theory, 193:105207, 2021.

[15]
F. Echenique, N. Immorlica, V.V. Vazirani, and A.E. Roth.
Online and Matching-Based Market Design.
Cambridge University Press, 2023.

[16]
Federico Echenique, Antonio Miralles, and Jun Zhang.
Constrained pseudo-market equilibrium.
American Economic Review, 111(11):3699–3732, November 2021.

[17]
Aris Filos-Ratsikas, Kristoffer Arnsfelt Hansen, Kasper Høgh, and Alexandros
Ppad-membership for problems with exact rational solutions: A general
approach via convex optimization, 2023.

[18]
Duncan K Foley.
Resource allocation and the public sector.(1967).

[19]
D. Gale and L. S. Shapley.
College admissions and the stability of marriage.
The American Mathematical Monthly, 69(1):9–15, 1962.

[20]
Yinghua He, Antonio Miralles, Marek Pycia, and Jianye Yan.
A pseudo-market approach to allocation with priorities.
American Economic Journal: Microeconomics, 10(3):272–314,

[21]
Mojtaba Hosseini and Vijay V. Vazirani.
Nash-bargaining-based models for matching markets: One-sided and
two-sided; fisher and arrow-debreu.
In Mark Braverman, editor, 13th Innovations in Theoretical
Computer Science Conference, ITCS 2022, January 31 - February 3, 2022,
Berkeley, CA, USA, volume 215 of LIPIcs, pages 86:1–86:20. Schloss
Dagstuhl - Leibniz-Zentrum für Informatik, 2022.

[22]
Aanund Hylland and Richard Zeckhauser.
The efficient allocation of individuals to positions.
Journal of Political economy, 87(2):293–314, 1979.

[23]
Nicole Immorlica, Brendan Lucier, Glen Weyl, and Joshua Mollner.
Approximate efficiency in matching markets.
In Nikhil R. Devanur and Pinyan Lu, editors, Web and Internet
Economics, pages 252–265, Cham, 2017. Springer International Publishing.

[24]
Jian Li, Yicheng Liu, Lingxiao Huang, and Pingzhong Tang.
Egalitarian pairwise kidney exchange: fast algorithms via linear
programming and parametric flow.
In Proceedings of the 2014 international conference on
Autonomous agents and multi-agent systems, pages 445–452, 2014.

[25]
Vikram Manjunath.
Fractional matching markets.
Games and Economic Behavior, 100:321–336, 2016.

[26]
Antonio Miralles and Marek Pycia.
Large vs. continuum assignment economies: Efficiency and
envy-freeness, 2016.

[27]
Ioannis Panageas, Thorben Tröbst, and Vijay V Vazirani.
Combinatorial algorithms for matching markets via nash bargaining:
One-sided, two-sided and non-bipartite.
arXiv preprint arXiv:2106.02024, 2021.

[28]
Alvin E Roth, Tayfun Sönmez, and M Utku Ünver.
Pairwise kidney exchange.
Journal of Economic theory, 125(2):151–188, 2005.

[29]
Hal Varian.
Equity, envy, and efficiency.
Journal of Economic Theory, 9(1):63–91, 1974.

[30]
Vijay V Vazirani.
The notion of a rational convex program, and an algorithm for the
arrow-debreu nash bargaining game.
Journal of the ACM (JACM), 59(2):1–36, 2012.

[31]
Vijay V. Vazirani and Mihalis Yannakakis.
Computational complexity of the Hylland-Zeckhauser mechanism for
one-sided matching markets.
In Innovations in Theoretical Computer Science, 2021.

[32]
John Von Neumann.
A certain zero-sum two-person game equivalent to the optimal
assignment problem.
Contributions to the Theory of Games, 2(0):5–12, 1953.

[33]
L. Zadeh.
Optimality and non-scalar-valued performance criteria.
IEEE Transactions on Automatic Control, 8(1):59–60, 1963.

[34]
Lin Zhou.
Strictly fair allocations in large exchange economies.
Journal of Economic Theory, 57(1):158–175, 1992.


# ALIVE

Matteo Cartiglia13, Arianna Rubino13, Shyam Narayanan1, Charlotte Frenkel 1,
Germain Haessig 2,
Giacomo Indiveri 1, Melika Payvand 1  1Institute of Neuroinformatics,
University of Zurich and ETH Zurich, Winterthurerstr. 190, 8057 Zurich,
Switzerland
2Austrian Institute of Technology Giefinggasse 4, 1210 Vienna, Austria 3 These
authors have contributed equally to this work.
Email: {rubinoa, camatteo<EMAIL_ADDRESS>
# Stochastic dendrites enable online learning in mixed-signal neuromorphic
processing systems

Matteo Cartiglia13, Arianna Rubino13, Shyam Narayanan1, Charlotte Frenkel 1,
Germain Haessig 2,
Giacomo Indiveri 1, Melika Payvand 1  1Institute of Neuroinformatics,
University of Zurich and ETH Zurich, Winterthurerstr. 190, 8057 Zurich,
Switzerland
2Austrian Institute of Technology Giefinggasse 4, 1210 Vienna, Austria 3 These
authors have contributed equally to this work.
Email: {rubinoa, camatteo<EMAIL_ADDRESS>
###### Abstract

The stringent memory and power constraints required in edge-computing sensory-
processing applications have made event-driven neuromorphic systems a
promising technology. On-chip online learning provides such systems the
ability to learn the statistics of the incoming data and to adapt to their
changes. Implementing online learning on event driven-neuromorphic systems
requires (i) a spike-based learning algorithm that calculates the weight
updates using only local information from streaming data, (ii) mapping these
weight updates onto limited bit precision memory and (iii) doing so in a
robust manner that does not lead to unnecessary updates as the system is
reaching its optimal output. Recent neuroscience studies have shown how
dendritic compartments of cortical neurons can solve these problems in
biological neural networks. Inspired by these studies we propose spike-based
learning circuits to implement stochastic dendritic online learning. The
circuits are embedded in a prototype spiking neural network fabricated using a
180 nm process. Following an algorithm-circuits co-design approach we present
circuits and behavioral simulation results that demonstrate the learning rule
features. We validate the proposed method using behavioral simulations of a
single-layer network with 4-bit precision weights applied to the MNIST
benchmark, and demonstrating results that reach accuracy levels above 85%.

###### Index Terms:

neuromorphic engineering, on-chip learning, dendritic processing, online
learning, HW-SW co-design

##  I Introduction

Our society is shifting to an era of pervasive specialized “edge computing”
for a wide variety of tasks. Artificial Intelligence (AI) is fueling this
trend by achieving remarkable results in pattern recognition. However, the
conventional computing technology used to run AI algorithms is not ideal for
edge-computing tasks that require minimal power consumption and always-on
real-time processing. Inspired by the computational principles of neural
processing systems, neuromorphic technologies have the potential to satisfy
the power-efficient and real-time processing requirements of edge-computing
applications [1]. The power efficiency of neuromorphic computing systems is
enabled by co-locating memory and processing [2] and by implementing fine-
grain parallelism strategies that remove the need to time-multiplex data
processing. This is achieved by instantiating multiple copies of the
processing elements (i.e., synapses and neurons) which have temporal dynamics
that are well-matched to those of the sensory signals of interest and which
operate in continuous physical time [3, 4]. Following this approach several
mixed-signal neuromorphic processing systems have been proposed that implement
spiking neural networks [5]. However, endowing these systems with online
learning abilities remains an open challenge.

Figure 1: Depiction of a pyramidal cortical neuron (left) and its multi-
compartment model used in artificial neural networks (right). Adapted from
[6].

Since providing a top-down error signal has been very successful in deep
learning[7, 8], some neuromorphic implementations have recently focused on
using errors for online learning [9, 10]. Moreover, many neuroscientific
studies have recently focused on error-based learning in the brain [11, 12].
In these studies, an error at the output units is assessed by comparing the
self-generated activity of the neurons with a target activity. This error is
then sent to the neurons as a feedback signal. Neurons are modeled with
multiple compartments with several active dendritic branches, each directly
linked to a somatic compartment (Fig. 1). The feedback signals arrive on the
distal Apical dendrites, while proximal Basal dendrites receive feed-forward
sensory information. Such types of models have been shown to implement error-
backpropagation using biologically-plausible plasticity mechanisms [11, 12,
6].

Figure 2:  (a) A layout-view of the 3x5 $mm^{2}$ design. The circuits proposed
are highlighted in the image and are composed of 4 individual neurons. (b) A
close-up of a single neuron. The dendritic branches, the soma, and the
learning block are all highlighted both in the schematic view and in the
physical layout. The Sensory input is made up of 40 excitatory and 16
inhibitory plastic synapses. The Teach input is made up of 3 excitatory and 5
inhibitory non-plastic synapses.

Inspired by these studies, we developed mixed-signal electronic circuits that
implement power-efficient online plasticity mechanisms that can be interfaced
to neuromorphic synapse and neuron circuits. Specifically, since the synaptic
weight precision is limited by the area of the memory that can be implemented
on-chip, an interface circuitry is required to translate the analog weight
updates to these low-bit precision memories. To solve this, we implement
stochastic rounding on chip, which converts the analog weight updates to a
probability of synaptic weight change [13, 14]. Additionally, this mechanism
provides robustness to the weight update process by reducing the probability
of updates as the error decreases to near-zero values.

We describe the software-hardware co-design approach that validated the
circuit designs with behavioral simulations, and that led to the fabrication
of the circuits using a standard 180 nm Complementary Metal-Oxide-
Semiconductor (CMOS) process (Fig. 2). In the next sections we present the
learning algorithm, the circuits implementations, and finally, the results
from the hardware and software simulations.

##  II The algorithm

###  II-A The multi-compartment neuron model

The selected neuron model used is a multi-compartment Leaky Integrate and Fire
(LIF) neuron in which currents from two dendritic branches, the bottom-up
Sensory branch ($I_{Sensory}$) and the top-down Teacher branch ($I_{Teach}$),
get integrated into a somatic compartment ($I_{Soma}$). This is summarized in
Eq. 1.

$\tau_{Soma}\cfrac{dI_{Soma}}{dt}+I_{Soma}(t)={{I_{Sensory}(t)}+{I_{Teach}(t)}}$
(1)

where $\tau_{Soma}$ is the neuron’s time constant.

The Sensory and Teach dendritric dynamics are modeled by Eq. 2 and 3
respectively.

$\tau_{Sensory}\cfrac{d{I_{Sensory_{j}}}}{dt}+{I_{Sensory_{j}}}=\sum\limits_{i}\sum\limits_{n}{I_{w_{ij}}}\delta_{i,s}(t-t_{n})$
(2)

$\tau_{Teach}\cfrac{d{I_{Teach_{j}}}}{dt}+{I_{Teach_{j}}}=\sum\limits_{n}{I_{w_{Teach}}}\delta_{i,t}(t-t_{n})$
(3)

where $i,j$ represents the input and output neurons respectively, and $n$ is
the time at which the input neurons spike. The $\delta_{i,s}$ and
$\delta_{i,t}$ function describes the input spike trains at the Sensory and
Teach branches, respectively. $\tau_{Sensory}$ and $\tau_{Teach}$ represent
the respective Sensory and Teach dendritic time constants, $I_{w_{ij}}$ is the
plastic synaptic strength between Sensory input $i$ and output neuron $j$
(Basal input), and ${I_{w_{Teach}}}$ is the non-plastic constant teacher
current (Apical input).

###  II-B The learning rule

The learning rule follows a dendritic-based prediction learning [15] which
resembles Delta rule. The Delta rule is the simplest form of gradient-descent
learning, minimizing a single-layer neural network’s Least Mean Square (LMS)
error. Here, the error is defined as the difference between the Sensory and
Teach dendritic branches as is described in Eq. 4.

${err_{ij}}=\alpha({I_{Teach_{j}}}-{I_{Sensory_{j}}})x_{i}|I_{Teach_{j}}|$ (4)

where $\alpha$ is the learning rate and the error is calculated every time a
pre-synaptic spike occurs ($x_{i}$). This error is scaled by the absolute
value of$(I_{Teach_{j}}$, which is in case of a true targets, is positive and
negative otherwise. In this way if the teaching signal is not present (i.e.,
$I_{Teach_{j}}=0$) the learning will stop.

To manage the changes on the limited-bit precision synaptic memory (here 4
bits) with the analog weight updates, that are proportional to the error
signal, we employed a stochastic rounding mechanism [13, 14, 16]. This is
implemented by comparing the magnitude of the error to a random number (here 6
bits), and if greater, incrementing or decrementing the weight depending on
the sign of the error. Otherwise, no weights are updated. Therefore, the
probability of weight update is proportional to the error in the network. This
is summarized in Eq. 5.

$\Delta w_{ij}(@x_{i})=\begin{cases}\pm
1\quad\text{if}\quad|err_{ij}|>rand()\\\ 0\quad\text{otherwise}\end{cases}$
(5)

Figure 3: Block diagram of the online learning mechanism. It consists of error
generation and stochastic rounding.

##  III The neuromorphic system

###  III-A The network architecture

The circuits have been fabricated in a prototype chip using a 180 nm CMOS
technology. The network occupies an area of $15mm^{2}$ ($5\,mm$ x $3mm$). The
total area occupied by the neural core presented here is $1.01mm^{2}$. The
physical layout of the network is shown in Fig. 2a.
The neural core can be broken down into four identical rows, where each row
has one neuron and 64 plastic and non-plastic synapses (Fig. 2b). The total
synaptic activity consists of two main dendritic branches. The Sensory branch
(Basal, ${I_{Sensory}}$) has 56 independent excitatory and inhibitory plastic
synapses (4 bit resolution), while the Teach branch (Apical, ${I_{Teach}}$),
has 8 excitatory and inhibitory non-plastic synapses. The excitatory synapses
inject a positive current into the dendrite, while the inhibitory ones takes
current away from it. Therefore, the total current in both dendritic branches
is the subtraction between the excitatory and inhibitory currents.

Upon the arrival of any input event on the Sensory dendritic branch, the Soma
integrates the weighted input, and if it reaches a threshold, it generates an
output spike. Moreover, the Sensory input event activates the learning block
which computes the error, generates a random number, and sends update signals
if the former is greater than the latter.

###  III-B The learning circuits

Figure 3 illustrates the block diagram of the learning mechanism, which
consists of the error-generation and stochastic rounding blocks.

####  III-B1 Error generation

The error generation implements Eq. 4, getting ${I_{Teach}}$ and
${I_{Sensory}}$ as inputs. The magnitude and the direction of weight update
are determined by the absolute value ($abs$), and Winner-Takes-All (WTA)
circuit, respectively.

The $abs$ circuit calculates the error by subtracting ${I_{Teach}}$ and
${I_{Sensory}}$ from each other in two directions, using $G_{1}$ and $G_{2}$
transconductance amplifiers, giving rise to the positive (${V_{PosErr}}$) and
negative (${V_{NegErr}}$) errors. Depending on which input is greater, either
the $M_{1}$ or the $M_{2}$ transistor is on, whose current flows in
$I_{AbsErr}$ giving rise to the absolute value of the error [17]. The $abs$
circuit is biased with ${I_{Teach}}$ to gate the generation of the error
signals only when the teacher is presented. If the teacher is not present
(${I_{Teach}}$ = 0), the $abs$ circuit is OFF. The Teach input can be either
positive or negative, however the $abs$ circuit is biased always with the
unsigned ${I_{Teach}}$. The learning block is designed to distinguish between
a positive and a negative Teach input and change the synaptic weights
accordingly.

The WTA decides the direction of update ($UP$ or $DOWN$) depending on which
one of ${V_{PosErr}}$ and ${V_{NegErr}}$ is higher [18]. If ${V_{PosErr}}$ and
${V_{NegErr}}$ are equal, $UP$ and $DOWN$ output of the WTA circuit are both
set to low, and thus no update signal is generated. This makes the learning
more robust by avoiding unnecessary updates as the system is reaching its
optimal output.

####  III-B2 Stochastic rounding circuits

The stochastic rounding block implements Eq. 5. Every time a pre-synaptic
spike is generated, a 6-bit pseudo-random number is generated from an on-chip
Linear Feedback Shift Register (LFSR) and is converted into an analog current
(${I_{random}}$) using a 6-bit current Digital to Analog Converter (DAC). The
${I_{random}}$ current is then compared to the error current, $I_{absolute}$,
generated by the error generation block (see III-B1). The comparison happens
through a modified version of Traff’s current comparator (CC) proposed in
[19]. Traff’s comparator input voltage is determined by the difference between
${I_{absolute}}$ and ${I_{random}}$ and is initially around ${Vdd/2}$. As soon
as ${I_{absolute}}$ and ${I_{random}}$ are not equal, the input voltage
increases or decreases and the change is amplified by the CC. The output
${OktoUpdate}$ is set to high/low if $I_{absolute}$ is larger/less than
$I_{random}$. The ${OktoUpdate}$ signal allows the weight update signals ($UP$
and $DOWN$) to pass through, only when it is high.

Since the probability of generating ${OktoUpdate}$ is proportional to the
error, as the error is reduced, the probability of generating an update is
also reduced. This leads to less weight updates as the system is reaching the
optimal output which makes the learning more robust and power efficient.

It is worth noting that the event-driven nature of the input allows us to use
a single stochastic rounding block for the entire row and generate an update
signal on the arrival of any pre-synaptic sensory input events.

##  IV Results

###  IV-A Circuit simulations

To validate the learning rule at a schematic-level, we show the simulation
results from one of the four identical rows of the chip (Figure 4). Initially,
a positive teacher (dotted blue in Fig. 4a) is presented, which causes the
weights of the excitatory synapses (dashdotted in orange in Fig. 4a) to
strengthen and minimize the error between the teacher and the sensory branch
(shown in red in Fig. 4a). At t = 0.3s, a negative teacher is shown to the
network, which causes the strengthened excitatory synapses to weaken, and the
inhibitory synapses (dashdotted in green in Fig. 4b) to strengthen and
minimize the error between the two dendritic branches.

Table I shows a direct comparison with previous works using on-chip on-line
error-based learning [16, 20]. The area in [16] has been normalized to a 180
nm technology node to allow for a better comparison with our work. The energy
per update in [20] is not directly comparable to our work, since the weight
update is calculated in a microcontroller, and the reported figure does not
take that into account. The 720 pJ reported in this work, is the total energy
in one row consumed to calculate and update one bit of synaptic memory.

Feature | [16] | [20] | This work
---|---|---|---
Resources | C5×5@10 | Re-configurable | 4 neurons
or topology | –FC128–FC10 |  | 64 synapse
Implementation | Digital | Digital | Mixed-signal
Learning | on-line | on-line | on-line
| error-based | re-configurable | error-based
Synaptic |  |  |
resolution | 8-bit | 1-9 bits | 4-bit
Static power | 61$\mu{W}$ | not-applicable11footnotemark: 1 | 455 nW
consumption |  |  | per neuron
Energy metric | N/A | 120 pJ11footnotemark: 1 | 720 pJ
|  | per weight update | per weight update
Technology node | 28 nm FDSOI | 14 nm FinFET | 180 nm CMOS
Power supply | 0.6 V | 0.75 V | 1.8 V
Area |  $13.2\,mm^{2}$ 22footnotemark: 2 | not-applicable 11footnotemark: 1 | $15\,mm^{2}$

[1]The learning block is implemented using a micro controller: area and power
figures are not directly comparable.
[2]Design normalized to a 180 nm technology node.

TABLE I: Comparison with previous works using on-chip on-line error-based
learning.

(a)

(b)

Figure 4: Cadence simulations of the learning signals of a single neuron; (a)
example of the signals involved in the learning process. The sensory dendritic
current (in red) is the sum of the excitatory (in orange) and inhibitory
currents (in green). A positive teacher signal (dotted in blue) is presented
until t = 0.3s. A negative teacher signal (dotted in purple) is presented
after that. (b) weight changes associated with the learning of (a). The
weights change their value to minimize the error between the sensory current
and the teacher current. An excitatory (dashdotted in orange) and an
inhibitory (dashdotted in green) weight are shown.

###  IV-B System

(a)

(b)

Figure 5: MNIST benchmarking results. (a) Accuracy comparison of a 784-10 SNN
network trained using surrogate gradients (SG) versus the same network trained
using the proposed stochastic dendritic online (SDO) rule. 80-20 randomized
division was used between training and testing. (b) Topology of the network in
the SDO simulations: input pixels are connected to the output layer with
excitatory (in orange) and inhibitory (in green) connections.

The proposed learning rule and system were fully modeled using the BRIAN2
simulator [21]. Fig. 5 shows the results of a shallow network (with topology
784-10 neurons) on the MNIST handwritten digit dataset. It shows a direct
comparison with a Spiking Neural Network (SNN) trained with surrogate gradient
descent [22] using 32-bit floating point (SG SNN), while the proposed
implementation uses 4-bit weights. The SG SNN reaches 91% accuracy while the
stochastic dendritic online implementation reaches above 85%. Interestingly
[23] achieves similar results with the same network topology.

Although both are SNNs, the difference in performance can be attributed to
that (i) the hardware network implements weight updates on streaming data,
with a batch size of one, compared to a batch size of 256 in SG SNN, and (ii)
the hardware uses only 4 bits of precision compared to the floating point
resolution in the SG SNN case.

##  V Discussion and Conclusion

We presented a prototype neuromorphic chip designed and fabricated in 180 nm
technology which implements stochastic dendritic-based online learning. We
proposed an algorithm–circuits co-design approach and validated it with
circuit simulations that demonstrate the proper operation of the learning
rule.

We benchmarked our system via system-level simulations, and obtained above 85%
accuracy on the MNIST dataset using only a one-layer network with 4-bit weight
precision.

As our circuits show that single neurons can adapt based on their incoming
activity, they represent a valid candidate for adaptive “edge computing”
applications that require online learning. To this end, future work will
include interfacing the chip with signals that need to be classified at the
edge, such as biomedical or industrial signals.

## Acknowledgment

The authors would like to thank Junren Chen for support during the last steps
in the chip layout. This project has received funding from the European
Union’s H2020 research and innovation programme under the H2020
BeFerrosynaptic (871737), MeM-Scales project (871371) and Sinergia Project
(CRSII5-180316).

## References

  * [1] G. Indiveri and S.-C. Liu.  Memory and information processing in neuromorphic systems.  Proceedings of the IEEE, 103(8):1379–1397, 2015.
  * [2] Abu Sebastian, Manuel Le Gallo, Riduan Khaddam-Aljameh, and Evangelos Eleftheriou.  Memory devices and applications for in-memory computing.  Nature nanotechnology, 15(7):529–544, 2020.
  * [3] Arianna Rubino, Can Livanelioglu, Ning Qiao, Melika Payvand, and Giacomo Indiveri.  Ultra-low-power FDSOI neural circuits for extreme-edge neuromorphic intelligence.  IEEE Transactions on Circuits and Systems I: Regular Papers, 68(1):45–56, 2020.
  * [4] Giacomo Indiveri and Yulia Sandamirskaya.  The importance of space and time for signal processing in neuromorphic agents.  IEEE Signal Processing Magazine, 36(6):16–28, 2019.
  * [5] Chetan Singh Thakur, Jamal Lottier Molin, Gert Cauwenberghs, Giacomo Indiveri, Kundan Kumar, Ning Qiao, Johannes Schemmel, Runchun Wang, Elisabetta Chicca, Jennifer Olson Hasler, Jae-sun Seo, Shimeng Yu, Yu Cao, André van Schaik, and Ralph Etienne-Cummings.  Large-scale neuromorphic spiking array processors: A quest to mimic the brain.  Frontiers in Neuroscience, 12:891, 2018.
  * [6] Jordan Guerguiev, Timothy P Lillicrap, and Blake A Richards.  Towards deep learning with segregated dendrites.  Elife, 6:e22901, 2017.
  * [7] Brian Crafton, Abhinav Parihar, Evan Gebhardt, and Arijit Raychowdhury.  Direct feedback alignment with sparse connections for local learning.  arXiv, pages 1–13, 2019.
  * [8] Friedemann Zenke and Emre O Neftci.  Brain-Inspired Learning on Neuromorphic Substrates.  Proceedings of the IEEE, pages 1–16, 2021.
  * [9] Melika Payvand, Mohammed E. Fouda, Fadi Kurdahi, Ahmed M. Eltawil, and Emre O. Neftci.  On-Chip Error-Triggered Learning of Multi-Layer Memristive Spiking Neural Networks.  IEEE Journal on Emerging and Selected Topics in Circuits and Systems, 10(4):522–535, 2020.
  * [10] Matteo Cartiglia, Germain Haessig, and Giacomo Indiveri.  An error-propagation spiking neural network compatible with neuromorphic processors.  In 2nd IEEE International Conference on Artificial Intelligence Circuits and Systems (AICAS), pages 84–88. IEEE, 2020.
  * [11] João Sacramento, Rui Ponte Costa, Yoshua Bengio, and Walter Senn.  Dendritic cortical microcircuits approximate the backpropagation algorithm.  In Advances in Neural Information Processing Systems, pages 8721–8732, 2018.
  * [12] James CR Whittington and Rafal Bogacz.  An approximation of the error backpropagation algorithm in a predictive coding network with local hebbian synaptic plasticity.  Neural computation, 29(5):1229–1262, 2017.
  * [13] Lorenz K Muller and Giacomo Indiveri.  Rounding methods for neural networks with low resolution synaptic weights.  arXiv preprint arXiv:1504.05767, 2015.
  * [14] Matthieu Courbariaux, Yoshua Bengio, and Jean-Pierre David.  Binaryconnect: Training deep neural networks with binary weights during propagations.  In Advances in neural information processing systems, pages 3123–3131, 2015.
  * [15] Robert Urbanczik and Walter Senn.  Learning by the dendritic prediction of somatic spiking.  Neuron, 81(3):521–528, 2014.
  * [16] Charlotte Frenkel, Jean-Didier Legat, and David Bol.  A 28-nm convolutional neuromorphic processor enabling online learning with spike-based retinas.  In 2020 IEEE International Symposium on Circuits and Systems (ISCAS), pages 1–5. IEEE, 2020.
  * [17] Melika Payvand, Manu V Nair, Lorenz K Müller, and Giacomo Indiveri.  A neuromorphic systems approach to in-memory computing with non-ideal memristive devices: From mitigation to exploitation.  Faraday Discussions, 213:487–510, 2019.
  * [18] M. Payvand and G. Indiveri.  Spike-based plasticity circuits for always-on on-line learning in neuromorphic systems.  In IEEE International Symposium on Circuits and Systems (ISCAS), pages 1–5. IEEE, 2019.
  * [19] H. Träff.  Novel approach to high speed cmos current comparators.  Electronics Letters, 28:310–311(1), January 1992.
  * [20] Mike Davies, Narayan Srinivasa, Tsung-Han Lin, Gautham Chinya, Yongqiang Cao, Sri Harsha Choday, Georgios Dimou, Prasad Joshi, Nabil Imam, Shweta Jain, Yuyun Liao, Chit-Kwan Lin, Andrew Lines, Ruokun Liu, Deepak Mathaikutty, Steven McCoy, Arnab Paul, Jonathan Tse, Guruguhanathan Venkataramanan, Yi-Hsin Weng, Andreas Wild, Yoonseok Yang, and Hong Wang.  Loihi: A neuromorphic manycore processor with on-chip learning.  IEEE Micro, 38(1):82–99, 2018.
  * [21] D. Goodman and R. Brette.  Brian: a simulator for spiking neural networks in Python.  Frontiers in Neuroinformatics, 2, 2008.
  * [22] Emre O Neftci, Hesham Mostafa, and Friedemann Zenke.  Surrogate gradient learning in spiking neural networks: Bringing the power of gradient-based optimization to spiking neural networks.  IEEE Signal Processing Magazine, 36(6):51–63, 2019.
  * [23] Charlotte Frenkel, Martin Lefebvre, Jean-Didier Legat, and David Bol.  A 0.086-mm2 12.7-pj/SOP 64k-synapse 256-neuron online-learning digital spiking neuromorphic processor in 28-nm CMOS.  IEEE Transactions on Biomedical Circuits and Systems, 13(1):145–158, 2019.

  *[CMOS]: Complementary Metal-Oxide-Semiconductor
  *[LIF]: Leaky Integrate and Fire
  *[LFSR]: Linear Feedback Shift Register
  *[DAC]: Digital to Analog Converter
  *[SNN]: Spiking Neural Network
  *[SNNs]: Spiking Neural Network


This paper has been withdrawn by the author due to the paper is far from
complishment.