# Free-Space Optical Spiking Neural Network

Reyhane Ahmadi
Department of Computer Engineering
Sharif University of Technology
Tehran, Iran, 11155-4365
<EMAIL_ADDRESS>
&Amirreza Ahmadnejad
Department of Electrical Engineering
Sharif University of Technology
Tehran, Iran, 11155-4365
<EMAIL_ADDRESS>
&Somayyeh Koohi
Department of Computer Engineering
Sharif University of Technology
Tehran, Iran, 11155-4365
<EMAIL_ADDRESS>

###### Abstract

Neuromorphic engineering has emerged as a promising avenue for developing
brain-inspired computational systems. However, conventional electronic AI-
based processors often encounter challenges related to processing speed and
thermal dissipation. As an alternative, optical implementations of such
processors have been proposed, capitalizing on the intrinsic information-
processing capabilities of light. Within the realm of optical neuromorphic
engineering, various optical neural networks (ONNs) have been explored. Among
these, Spiking Neural Networks (SNNs) have exhibited notable success in
emulating the computational principles of the human brain. Nevertheless, the
integration of optical SNN processors has presented formidable obstacles,
mainly when dealing with the computational demands of large datasets. In
response to these challenges, we introduce a pioneering concept: the Free-
space Optical deep Spiking Convolutional Neural Network (OSCNN). This novel
approach draws inspiration from computational models of the human eye. We have
meticulously designed various optical components within the OSCNN to tackle
object detection tasks across prominent benchmark datasets, including MNIST,
ETH 80, and Caltech. Our results demonstrate promising performance with
minimal latency and power consumption compared to their electronic ONN
counterparts. Additionally, we conducted several pertinent simulations, such
as optical intensity-to-latency conversion and synchronization. Of particular
significance is the evaluation of the feature extraction layer, employing a
Gabor filter bank, which stands to impact the practical deployment of diverse
ONN architectures significantly.

##  1 Introduction

The human brain represents a profoundly intricate and remarkable biological
entity. The endeavor to engineer a computational processor possessing
commensurate attributes in power, precision, integration, and speed has
perennially constituted a paramount aspiration for processor designers.
Neuromorphic Engineering (NE) [1] stands as a foundational paradigm
facilitating the realization of such processors, primarily through the
incorporation of neural network architectures (NNs). Despite the notable
achievements resulting from this approach [2, 3], the central challenge in
processor design endures as the demand for processing voluminous datasets
continues to burgeon. To address this persistent challenge, optical
Neuromorphic Engineering has emerged as a novel and innovative domain. Optical
Neuromorphic Engineering exploits the distinctive attributes of light, which
include its exceptional propagation speed and the extended degrees of freedom
it affords in comparison to electrons, encompassing characteristics like path,
frequency, phase, polarization, and mode. Furthermore, optical systems
manifest reduced loss, rendering them remarkably compelling for the design and
construction of Optical Neural Networks (ONNs).

The pursuit of processor miniaturization has perennially remained a core
objective in the field of processor design [4]. This pursuit has extended to
the domain of Optical Neural Networks (ONNs) [5]. Nevertheless, the efficient
processing of substantial datasets at elevated speeds presents a considerable
challenge within this sphere. In response to this challenge, recent research
endeavors have refocused on developing processors harnessing the capabilities
of optical free space (OFS) devices [6, 7, 8, 9, 10].

Spiking Neural Networks (SNNs), constituting a class of neural networks that
emulate the structural and functional aspects of the human brain, have
garnered significant attention in this context. Many optical models have been
proposed for implementing SNNs; however, these models have been primarily
integrated. Several intricate designs featuring components such as Vertical-
Cavity Surface-Emitting Lasers (VCSELs) [11], micro ring resonators [12], and
phase-change materials [13] have been suggested. Nonetheless, these designs
prove ill-suited for the demanding task of high-volume data processing,
thereby underscoring the formidable challenges encountered within the realm of
optical Neuromorphic Engineering.

To develop an Optical Free Space (OFS) model for Spiking Neural Networks
(SNNs), it is imperative to establish precise mathematical models for each
constituent component. Remarkably successful models rooted in neuroscience
[14, 15] have been devised, emulating the structural attributes responsible
for object detection within the human eye. We aim to draw upon these well-
established models as a source of inspiration for designing OFS components,
thereby facilitating the simulation of the Free-Space Optical deep Spiking
Convolutional Neural Network (OSCNN). To the best of our knowledge, OSCNN
marks the inaugural foray into the realm of OFS modeling for SNNs,
encompassing critical elements such as Gabor filters for feature extraction,
intensity-to-delay conversion, synchronization mechanisms, convolution layers,
Max-pooling procedures, and a classification framework.

To execute the intensity-to-delay conversion, a dedicated module was
introduced, employing a Spatial Light Modulator (SLM) after the feature
extractor layer. Moreover, an optical synchronizer has been meticulously
devised to address the temporal processing aspects inherent to optic signals.
Ensuring that the time order of signals remains intact post-convolution, this
synchronizer draws inspiration from the Free-Space Optical (FSO) delay line
concept [16].

The performance of OSCNN was systematically evaluated across three distinct
datasets: MNIST, Caltech, and ETH80. OSCNN demonstrated notable achievements,
boasting significant performance metrics compared to electronic Neural
Networks (NNs) and Optical Neural Networks (ONNs). Notably, Gabor filters were
harnessed as feature extractors in the initial model layer, with evaluations
conducted under both trained and fixed conditions. While alternative filters,
such as Canny, Laplacian, and Sobel, were explored as feature extraction
mechanisms, the most favorable outcomes for OSCNN were attained using Gabor
filters. The results underscore the versatility of Gabor-form convolutional
kernels, revealing their efficacy in image and time-series processing
applications [17, 18].

##  2 OSCNN Model

In delineating the architectural framework of the Optical Free Space Spiking
Convolutional Neural Network (OSCNN), it is imperative to establish a
comprehensive model for neurons within Spiking Neural Networks (SNNs, also
known as integrate-and-fire models). SNNs fundamentally operate on the
principles of spike-timing-dependent plasticity (STDP), a paradigm
necessitating optical modeling to faithfully replicate its mechanisms. OSCNN
adopts a specialized variation of computational neurons, acknowledging that
forthcoming research will delve into developing a more precise neuron model.
It is worth noting that as an alternative to STDP, backpropagation (BP) can be
employed for training the OSCNN model, as described in [19]. Subsequently,
this section elaborates on the mathematical underpinnings of each module
within OSCNN, along with their optical equivalents. The holistic structure of
the OSCNN model is illustrated in Figure 1 for a comprehensive overview.

Figure 1: Outline of OSCNN.

###  2.1 Gabor Filters

A Gabor filter serves as a bandpass filter meticulously designed to emulate
the receptive field characteristics of neurons in the primary visual cortex of
the human brain [20]. This filter is uniquely characterized by two key
parameters: the frequency of the sinusoidal wave and the width of the Gaussian
envelope. In implementing feature extraction via convolution using the Gabor
kernel filter, an optical 4f correlator emerges as a viable solution [7, 10].

To execute the convolution of an image, denoted as $x$, with a Gabor filter
represented by $g$, a Fourier domain approach is adopted. In this process, the
Gabor filter and the input image are transformed into their respective Fourier
domain representations. Subsequently, these transformed representations are
multiplicatively combined within the Fourier domain. The Fourier transform of
the Gabor filter is expressed as follows:

$\begin{split}G(u,v)=e^{-\frac{1}{2}\left(\frac{\sigma_{x}^{2}(u-k_{x})^{2}+\sigma_{y}^{2}(v-k_{y})^{2}}{\sigma_{x}^{2}\sigma_{y}^{2}}\right)}\end{split}$
(1)

The Fourier transformation of the Gabor filter reveals itself as another Gabor
function, characterized by specific frequency components, denoted as $k_{x}$
and $k_{y}$, alongside standard deviations represented by $\sigma_{x}$ and
$\sigma_{y}$. Consequently, the mathematical model for the convolution of an
image $X(u,v)$ with a Gabor filter through the utilization of an optical 4f
correlator can be expressed in the following manner:

$\begin{split}y_{1}(u,v)=X(u,v)\cdot
G(u,v)){\sigma_{x}^{2}\sigma_{y}^{2}}\end{split}$ (2)

###  2.2 Intensity to Phase Conversion

The resultant features manifest as light patterns characterized by varying
intensities. To faithfully represent the optical properties and mitigate
signal losses, a crucial transformation process is invoked, converting the
information modulated on intensity in the input image into phase information.
This transformation is known as intensity-to-phase (latency) conversion and is
effectuated through a Spatial Light Modulator (SLM).

Suppose $y_{1}(u,v)$ symbolizes the output derived from the preceding
convolution module, associated with intensity values $I(u,v)$. In that case,
the conversion operation can be articulated as
$I(u,v)=I_{0}\cos^{2}(\phi(u,v))$, where $I_{0}$ designates the maximum
intensity and $\phi(u,v)$ signifies the phase values. To effectuate this
conversion, the SLM introduces a spatially varying phase shift commensurate
with the intensity at each pixel. Mathematically, the SLM can be represented
as a complex transmission function $T(u,v)$, satisfying the condition
$|T(u,v)|\leq 1$. In most instances, $T(u,v)$ is expressed as
$e^{i\phi(u,v)}$, where $\phi(u,v)$ corresponds to the desired phase shift at
each pixel. Hence, the complex transmission function assumes the form of a
complex exponential, featuring a phase term directly proportional to the
desired phase shift. Consequently, the output stemming from the intensity-to-
phase conversion module is succinctly represented as:

$\begin{split}{y}_{2}(u,v)=y_{1}(u,v)\cdot T(u,v)\end{split}$ (3)

###  2.3 Optical Synchronizer

Upon converting intensity to phase, it is imperative to address that each
optical signal propagates at varying speeds, with the fastest signal being of
utmost significance. To facilitate the concurrent processing of all signs and
their combinations, it is essential to ensure their emission commences from a
common temporal reference point. This task is conventionally achieved in
integrated ONNs by employing a delay line structure [2]. In contrast, this
practice is less prevalent within Optical Free Space (OFS). In OFS, the
synchronization is typically accomplished through a diffraction grating and a
Spatial Light Modulator (SLM) instead of the traditional delay line or
parallel mirrors [9].

A synchronizer is engineered by leveraging an SLM and a diffraction grating by
the principles delineated in [16]. The output field originating from the
preceding module is $y_{2}(u,v)$, subjected to diffraction through a grating
characterized by a specific pitch denoted as $\Lambda$, along with a
designated diffraction angle, $\theta$. The output emerging from the grind
embodies multiple diffraction orders, each encompassing a version of the input
signal delayed by a distinct temporal offset. The electric field associated
with the $m$th diffraction order is aptly described as $E_{m}(u,v)$, featuring
a time delay of $\tau_{m}=m\Lambda\sin\theta/c$. Consequently, the electric
field affiliated with the $m$th diffraction order can be succinctly
articulated as:

$\begin{split}E_{m}(u,v)=y_{2}(u,v)e^{i2\pi m\Lambda sin\theta/c}\end{split}$
(4)

After the diffraction grating, the individual diffraction orders undergo phase
modulation through the utilization of an SLM, where the relative phases of the
input signals are meticulously adjusted. The phase modulation introduced by
the SLM can be mathematically represented by a complex-valued function
$f(u,v)$. To reconstitute the electric fields associated with the diverse
diffraction orders, a lens boasting a focal length denoted as $f$ is employed.
Ultimately, the output signal is derived from the electric field situated at
the focal point of the lens. This output signal can be expressed as (with
$\mathrm{sinc}(x)=\frac{\sin(\pi x)}{\pi x}$):

$\begin{split}y_{i}(u,v)=\sum_{m=-\infty}^{\infty}E_{m}(u,v)e^{if(u,v)}sinc(\frac{m\Lambda
sin\theta}{f})\end{split}$ (5)

###  2.4 Layers and Classifier

Following the feature extraction phase, the optical signals must be
amalgamated to capture the salient features embedded within the image. This
synthesis is achieved via a 3-layer system comprising a convolution layer, a
synchronizer, and a max-pooling module [15]. The convolution layer effectively
combines various optical signals, each subject to different learnable weights,
while the max-pooling layer identifies and preserves the most crucial
features. The necessity for a synchronizer arises from the typical propensity
of the convolution operation to disarrange the temporal order of signals. The
convolution operation is mathematically represented as a 4f correlator,
characterized by trainable kernels [7, 10]. In parallel, the max-pooling
function is realized by employing a 4f correlator, incorporating a saturable
absorber (SA) as a nonlinearity, which can be conceptually likened to a
multilayer neural network augmented with nonlinearity [7, 21].

The optical emulation of a classifier is inspired by the approach delineated
in [7]. It consists of a single classifier utilizing an MNN with a non-
linearity module, specifically the saturable absorber (SA). The comprehensive
structural layout of the OSCNN model is prominently depicted in Figure 2.

Figure 2: The simulation of OSCNN involves loading data from a laser and SLM
and simulating different optical free space components.

##  3 Simulation Study

This section delves into the comprehensive simulations conducted to evaluate
the performance of the Free-Space Optical Spiking Convolutional Neural Network
(OSCNN). The study extends to comparative analyses with other models,
encompassing both electrical and optical domains, encompassing free-space and
integrated approaches. A primary focus of the investigation lies in the in-
depth analysis of the first-layer feature extractor kernels in both fixed and
trainable configurations. Incorporating Gabor filters as convolutional kernels
in the feature extractors of Convolutional Neural Networks (CNNs) is
particularly emphasized, owing to its biological inspiration from the human
brain and inherent properties. This is further juxtaposed with comparisons
involving other well-recognized filters such as Sobel, Canny, and Laplacian.
Furthermore, the impact of noise on the input image is methodically examined,
and the temporal consumption of electrical or optical resources by the OSCNN
is quantitatively measured.

The OSCNN model was subjected to rigorous training using the MNIST dataset,
facilitated by a V100 Tesla GPU on the Google Colab platform. The cumulative
duration for processing and training with the MNIST dataset amounted to 2
hours and 37 minutes, a timeframe that is markedly consistent with processing
times associated with other electrical and optical models, as elucidated in
[7, 8, 22, 10]. The training process was executed through backpropagation,
functionally equivalent to the Spike-Timing-Dependent Plasticity (STDP)
process described in [19].

In addition to the MNIST dataset, the model was subjected to rigorous training
and testing on the ETH-80 dataset [23] and the Caltech dataset [24], both of
which are well-recognized benchmarks within the Spiking Neural Network (SNN)
domain [15]. The simulation of the OSCNN closely aligns with established
methodologies applied in prior models, such as those documented in [7, 8, 9,
25, 10]. Each optical module is rigorously formulated mathematically within
these simulations and is referred to as a behavioral model.

###  3.1 Gabor Filters

The initial layer of the OSCNN is characterized by an array of Gabor filters,
each possessing distinct spatial orientations and thicknesses. To provide a
comprehensive comparative analysis of the feature extraction process, we
benchmark the OSCNN against other Optical Neural Networks (ONNs). For
instance, the Diffractive Deep Neural Network (D2NN) [26] leverages light
diffraction properties employing apertures designed through the Huygens
principle. In the OSCNN, Gabor filters are deployed to extract the most
salient features embedded within the images meticulously. To this end, various
feature extraction approaches are examined, including fixed Gabor filters
devoid of training, trainable Gabor filters, and established filters like
Canny, Laplacian, and Sobel. These diverse feature extractors are employed
with the MNIST, Caltech, and ETH-80 datasets. The ensuing influence of these
filters on the output accuracy is meticulously assessed and is presented in
Table 1.

Table 1: OSCNN accuracy with different feature extractors Filter | MNIST | Caltech | ETH-80
---|---|---|---
Fixed Gabor | $91.4$ | $88.3$ | $84.7$
Trainable Gabor | $\bf{95.2}$ | $\bf{91.3}$ | $\bf{89.7}$
Canny | $89.2$ | $86.7$ | $84.2$
Laplacian | $90.4$ | $84.3$ | $83.5$
Sobel | $86.7$ | $82.1$ | $79.4$

The outcomes, as depicted in Table 1, reveal that the trainable Gabor filter,
endowed with adaptable parameters about its filter length and central
frequency, attains the highest performance among the filters tested. However,
it is noteworthy that even the fixed Gabor filter consistently outperforms the
alternative filters. This substantiates our assertion that Gabor filters can
be regarded as reliable and effective feature extractors across various
Optical Neural Networks (ONNs) and can be effectively deployed in the first
layer as Convolutional Neural Network (CNN) kernels. For further visual
representation, Figure 3 showcases the output images generated by applying the
Gabor input filter to an image from the dataset, exemplified by image number
8.

Figure 3: Gabor filters are created as sine wave filters in four different
directions, namely $0$, $\frac{\pi}{2}$, $\pi$, and $\frac{3\pi}{4}$ degrees.
These filters have different thicknesses and function as edge detectors,
resembling the simple cells found in the primary visual cortex [20].

This comprehensive analysis underscores the exceptional efficacy of Gabor
filters as optimal models for feature extraction in the initial layer of
Optical Neural Networks (ONNs). Consequently, Gabor filters stand as a highly
recommended choice for the processing of diverse data types, extending their
applicability to various domains, including image analysis and the handling of
temporal signals, such as audio, as expounded in [17].

###  3.2 Intensity to Phase Conversion

Employing a Spatial Light Modulator (SLM) for intensity-to-phase conversion
within OSCNN necessitates meticulous considerations to ensure optimal
performance. These considerations encompass spatial resolution, phase
modulation range, and the selection of the optimal operating wavelength. The
spatial resolution of the SLM plays a pivotal role in defining the output
spatial resolution of the OSCNN. Therefore, a high-resolution SLM is highly
desirable for achieving superior spatial resolution. Furthermore, the phase
modulation range exhibited by the SLM profoundly influences the dynamic range
of the OSCNN, with a more extensive phase modulation range facilitating a
broader dynamic range. The SLM’s specifications must align optimally with the
operating wavelength of the optical system employed within the OSCNN to ensure
maximum efficiency and precision.

While the simulations affirm that the SLM is one of the most straightforward
and practical choices for implementing Intensity Latency conversion,
alternative methods are available, including Amplitude-Phase Modulation (APM)
and holographic filtering. The selection among these methods should be
contingent on the specific requirements of the application, encompassing
factors such as speed and power consumption. To provide a comparative
assessment of processing speed for a $28\times 28$ pixel image with Gabor
filters, followed by intensity to latency conversion, the time required for
the SLM stands at approximately $1ms$ [27]. In contrast, for APM, it is
roughly $0.5ms$ [28], and for holographic filtering, it approximates 1 ms
[29]. Therefore, the speed comparison can be briefly summarized as follows:

$t_{APM}<t_{SLM}<t_{HF}$ (6)

Where $t_{HF}$ is the time required for intensity to latency conversion using
holographic filtering.

In power consumption, it is essential to note that SLMs primarily leverage
electro-optic effects, which typically demand power in the millijoules per
second range [30]. In contrast, Amplitude-Phase Modulation (APM) leans on
nonlinear effects like Pockels or Kerr, which may escalate power consumption
to the level of joules per second [28]. Therefore, the comparison of power
consumption can be concisely summarized as follows:

$P_{SLM}=P_{HF}<P_{AMP}$ (7)

###  3.3 Synchronizer

The accurate simulation of a synchronizer hinges on a meticulous planning
process, necessitating careful consideration of the diffraction angles and
grating period. This planning is crucial to ensure that the diffraction orders
are distinctly separated and that the resulting synchronized output field
boasts a superior signal-to-noise ratio. The Spatial Light Modulator (SLM)
plays a pivotal role in this endeavor by encoding the filter as a phase
pattern. Precise design of the phase pattern on the SLM, in conjunction with
careful spacing of the grating, enables the synchronization of distinct
diffraction orders, permitting their simultaneous arrival and concurrent
processing.

A diffraction grating, in this context, assumes the form of a periodic
structure replete with regularly spaced lines or slits, adept at diffracting
incident light into a specific pattern. When designing a diffraction grating,
two paramount parameters must be meticulously considered: the grating period,
denoted as $d$, and the number of lines, represented by $N$. To craft a
grating tailored for a specific wavelength, the following formula proves
instrumental:

$dsin(\theta)=m\lambda$ (8)

In our simulation of the OSCNN, we employ a red laser with a wavelength of
$1550nm$, and the separation between the SLM and the diffraction grating is
$5cm$. For this scenario, we consider diffraction of the incident light at an
angle of $5^{\circ}$ while specifically utilizing the first diffraction order,
$m=1$. As such, referring to equations 8 , we calculate the grating period to
be $11.8\mu m$, and the corresponding diffraction angle is approximately
$32.8^{\circ}$.

###  3.4 Model Performance

The evaluation of OSCNN’s proficiency in object detection tasks encompassed
its rigorous training and testing on three distinct datasets: MNIST, ETH 80,
and Caltech. Subsequently, the results were meticulously juxtaposed against
established electronic and optical models, both in the free-space and
integrated domains. These comparative results are illustrated in Figure 4 and
tabulated in Table 2.

Figure 4 offers a comparative overview, contrasting the performance of
electronic implementations, specifically those of Alexnet and electrical
models, with their OFS counterparts. It is noteworthy that SNNs, whether
instantiated electronically or optically, consistently manifest superior
object detection performance when compared to CNNs.

Figure 4: Accuracy of different electrical and optical Neural Network on Caltech and ETH-80 data sets. Table 2: MNIST Object Detection Accuracy and Models Implementation | Model | MNIST Accuracy (%) | Input Size
---|---|---|---
Electrical | Perception | MLP [31] | 92.3 | $28\times 28$
Convolution | AlexNet [22] | 98.3 | $28\times 28$
Spiking | Kheradpisheh et al. [15] | 97.2 | $28\times 28$
Optical | Perception | Shen et al. [3] | 95.4 | $3\times 3$
D2NN [23] | 99.1 | $3\times 3$
Ryou et al. [6] | 86.7 | $3\times 3$
Convolution | Bagherian et al. [2] | 87.7 | $3\times 3$
Sadeghzadeh et al. [7] | 97.6 | $28\times 28$
Spiking | Feldmann et al. [13] | 98.3 | $3\times 3$
| Xiang et al. [12] | 89.9 | $3\times 3$
| OSCNN | 95.2 | $28\times 28$

Table 2 serves as a comprehensive comparative analysis, shedding light on the
performance of various electrical and optical models, all tested on the MNIST
dataset. The findings unequivocally underscore that OSCNN exhibits commendable
performance levels in the context of object detection tasks. Nevertheless, it
is crucial to acknowledge a notable disparity between the implementations of
NNs and ONNs, a distinction prominently illustrated in Figure 4 and elaborated
upon in Table 2. In this context, it is essential to highlight that FSO
implementations consistently demonstrate accelerated processing speeds when
juxtaposed against their integrated counterparts. Furthermore, FSO
implementations exhibit remarkable data-handling capabilities, extending to
vast datasets, including biological data, a feat that remains presently
unattainable for integrated ONNs.

###  3.5 Noise Robustness

The assessment of OSCNN extends to its resilience against various levels of
input image noise. In our evaluation, we introduced white noise to the input
images, thereby subjecting the model to different noise levels ranging from
$5\%$ to $50\%$. Our findings regarding OSCNN’s performance under varying
noise conditions are vividly depicted in Figure 5. Furthermore, we subjected
several other electrical and optical Neural Networks (NNs) to evaluations to
gauge their robustness against noise. It is worth noting that up to a noise
level of $15\%$, the model demonstrates a reasonable degree of accuracy.
However, as noise levels surge beyond this threshold, accuracy experiences a
steep decline, ultimately converging to the chance level of $40\%$.

Figure 5: Accuracy of different electrical and optical Neural Networks for
various input noise levels.

Figure 4 offers a compelling visual representation of the comparative
robustness exhibited by electrical CNNs and Optical formats, notably
surpassing the robustness demonstrated by electrical SNNs and Optical ones.
This observation prominently underscores one of the primary limitations
inherent to SNNs, which is discernible not only in the electrical domain but
also extends into the optical domain. Consequently, it becomes evident that
further research endeavors are warranted to soothe and bolster the performance
of SNNs, particularly in the context of robustness against noise.

##  4 Speed Analysis

Optical computing inherently harbors several compelling capabilities. The
paramount advantage inherent to this optical implementation is its exceptional
capacity for accelerated processing, distinctly outpacing its electrical
counterparts. To provide a comprehensive insight into the relative speed of
OSCNN, an intricate assessment of latency can be explained through the
following estimation:

$\displaystyle t_{latency}$
$\displaystyle=t_{source}+t_{Gaborfilters}+t_{PM}+t_{Sync}$ (9)
$\displaystyle+t_{Conv1}+t_{Sync}+t_{Maxpooling}$
$\displaystyle+t_{Conv2}+t_{Sync}+t_{Maxpooling}$
$\displaystyle+t_{Conv3}+t_{Sync}+t_{Maxpooling}$
$\displaystyle+t_{classifier}+t_{camera}+t_{transferdata}$

In this comprehensive analysis, the estimation of OSCNN’s latency is
meticulously deconstructed, considering various pertinent factors. To begin,
$t_{source}$ signifies the modulation delay associated with input images. With
the consideration of Spatial Light Modulators (SLMs) featuring a 1 kHz
switching frequency, $t_{source}$ is approximated at $1ms$ [27]. Moving
forward, when addressing $t_{Gaborfilters}$, which fundamentally represents a
convolutional layer employing Gabor kernels, the estimated processing time is
approximately $5ps$, a minuscule interval that can be safely disregarded. As
for converting intensity into latency, a process executed by $t_{PM}$, the
approximate duration allocated to this operation is also pegged at around
$1ms$ [27].

The synchronization system encompasses both a grating layer and an SLM. The
dissonant layer’s processing speed is approximately $1ms$ for a $28\times 28$
image [16], and the cumulative time required for processing the
synchronization block is approximately $2ms$. Subsequently, we arrive at
$t_{Conv1}$ and $t_{Maxpooling}$, representing the optical propagation delays
within the convolution and pooling layers. When employing 4f optical
correlators, the processing delay for each 4f layer is approximately $5ps$,
equivalently, $5ps$ for the pooling layers [10], durations that can be
conveniently deemed negligible.

The classification module, integral to the system, comprises a 4F multiplier
and a nonlinear Saturable Absorber (SA) unit. The delay attributed to the
nonlinear unit is approximately $25ns$ [10], while the classification
component exhibits a delay of roughly $25ms$. In this context, $t_{camera}$
represents the time required for photodetectors to capture and convert output
images into electrical data. A high-speed commercial camera, capturing images
at a rate of 2500 frames per second [7], is characterized by an estimated
latency of $0.4ms$. Lastly, the communication interface introduces delays when
transmitting the camera’s output data to a computer, as reflected by
$t_{transferdata}$. By leveraging USB 3.1 Gen2 at a data rate of 10 Gbit/s and
processing a 50 kB image, $t_{transferdata}$ is quantified at $0.04ms$.
Consequently, the cumulative delay attributed to OSCNN is reasonably estimated
at approximately $2.44ms$.

##  5 Power Consumption Analysis

The research outlined in [10] sheds valuable light on the power dynamics
within the OSCNN architecture. It delineates that Convolution, nonlinearity,
and pooling operations are characterized by minimal energy consumption.
Consequently, it follows that the principal contributor to energy consumption
within the system lies in the domain of signal transmission [10]. To
contextualize these findings, it’s essential to consider that each pixel’s
capture entails a power demand of 1W, as noted in [7]:

$P_{optical}=\frac{n^{2}\times n_{kernel}}{\eta\times t^{p}}$ (10)

The unit of power in this context is expressed in micro-watts. Within the
equation, $n^{2}$ denotes the total number of pixels per 4f correlator system,
while $p$ signifies the number of optical elements traversing the optical
path. The variable $n_{kernel}$ is an arbitrary value denoting the number of
kernels utilized by the convolutional layer. The parameter $t$ represents the
fraction of incident power received by each optical element. Additionally,
$\eta$ embodies the source efficiency. Conversely, when endeavoring to
calculate the electric power consumption, careful consideration must be
extended to the following factors, as expounded in [10]:

$P_{electrical}=\beta\times n^{2}\times k^{2}\times n_{kernel}\times
P_{switch}$ (11)

In the equation presented, the variable $\beta$ is determined by the
architectural characteristics of the program, with $k$ signifying the kernel
size and $P_{switch}$ representing the energy consumed by each operation. It’s
imperative to acknowledge that as the kernel size expands, electronic
components invariably escalate their power consumption. Conversely, the
optical implementation of convolutional layers showcases a distinct trend: as
the kernel size diminishes, the power consumption concurrently decreases. This
divergence in behavior highlights a remarkable aspect of optical
implementations—namely, their capacity for significant power savings,
particularly in scenarios featuring large kernel sizes, in contrast to their
electrical counterparts.

It’s also instructive to juxtapose these power considerations with the energy
demands of the human brain during image processing. As evidenced by research
in [28], the aggregate power consumed by the human brain for image processing
hovers around 0.2 watts, accomplished within a processing timeframe of
approximately ten milliseconds. In light of these findings, it becomes
apparent that the optical neuromorphic approach, encompassing the optical
creation of DNNs and SNNs, holds the potential to propel us a significant step
closer to the realization of a processor endowed with computational prowess
akin to the human brain, marked by its remarkable power efficiency.

##  6 Conclusion

This article marks a significant milestone by introducing the pioneering
Optical Deep Spiking Convolutional Neural Network model operating in free
space. Drawing inspiration from the computational model of the human eye, this
model excels in detecting patterns with commendable accuracy, processing
speed, and power efficiency. An essential revelation from this study is the
suggestion that, in the realm of optical neural networks tasked with image
processing, the initial layer can be effectively realized by deploying Gabor
filters, thereby revolutionizing the approach to feature extraction. Moreover,
the remaining optical free-space components, encompassing the Intensity-to-
Delay conversion and the Synchronizer, were meticulously designed by
harnessing the potential of readily available optical components. It’s
important to note that this optical model refrains from delving into the
intricacies of biological neuron modeling and Spike-Timing-Dependent
Plasticity (STDP) training. These aspects are deliberately deferred to future
research endeavors, where the focus would entail designing a dedicated
structure based on resonators or topological photonics. The objective is to
authentically simulate the precise neuron model and subsequently introduce an
array structure to usher in novel designs for convolutional layers, akin to
the advancements witnessed in metasurface technology. Additionally, the study
relies on a well-established behavioral model instead of a numerical analysis
of electromagnetic fields. By introducing optical free space Deep Spiking
Convolutional Neural Network models, a significant stride is taken toward the
realization of high-powered, high-speed processors inspired by the human
brain. This endeavor propels us closer to the development of an artificial
brain, manifesting itself in an optical form endowed with formidable
computational capabilities.

##  7 Acknowledgments

We acknowledge support from Dr. Zahra Kavehvash (Electrical Engineering
Department of Sharif University of Technology) and Dr. Alireza Ejlali
(Computer Engineering Department of Sharif University of Technology). We also
would like to express our gratitude to OpenAI for the invaluable assistance
provided by ChatGPT [32] during the revision and refinement of this paper.

## References

  * [1] Paul R Prucnal and Bhavin J Shastri.  Neuromorphic photonics.  CRC press, 2017.
  * [2] Hengameh Bagherian, Scott Skirlo, Yichen Shen, Huaiyu Meng, Vladimir Ceperic, and Marin Soljacic.  On-chip optical convolutional neural networks.  arXiv preprint arXiv:1808.03303, 2018.
  * [3] Yichen Shen, Nicholas C Harris, Scott Skirlo, Mihika Prabhu, Tom Baehr-Jones, Michael Hochberg, Xin Sun, Shijie Zhao, Hugo Larochelle, Dirk Englund, et al.  Deep learning with coherent nanophotonic circuits.  Nature photonics, 11(7):441–446, 2017.
  * [4] Robert W Keyes.  Miniaturization of electronics and its limits.  IBM Journal of Research and Development, 32(1):84–88, 1988.
  * [5] Runqin Xu, Pin Lv, Fanjiang Xu, and Yishi Shi.  A survey of approaches for implementing optical neural networks.  Optics & Laser Technology, 136:106787, 2021.
  * [6] Albert Ryou, James Whitehead, Maksym Zhelyeznyakov, Paul Anderson, Cem Keskin, Michal Bajcsy, and Arka Majumdar.  Free-space optical neural network based on thermal atomic nonlinearity.  Photonics Research, 9(4):B128–B134, 2021.
  * [7] Hoda Sadeghzadeh and Somayyeh Koohi.  High-speed multi-layer convolutional neural network based on free-space optics.  IEEE Photonics Journal, 14(4):1–12, 2022.
  * [8] Hoda Sadeghzadeh and Somayyeh Koohi.  Translation-invariant optical neural network for image classification.  Scientific Reports, 12(1):17232, 2022.
  * [9] Aida Ebrahimi Dehghanpour and Somayyeh Koohi.  All-optical recurrent neural network with reconfigurable activation function.  IEEE Journal of Selected Topics in Quantum Electronics, 29(2: Optical Computing):1–14, 2022.
  * [10] Hoda Sadeghzadeh, Somayyeh Koohi, and Ali Fele Paranj.  Free-space optical neural network based on optical nonlinearity and pooling operations.  IEEE Access, 9:146533–146549, 2021.
  * [11] Shuiying Xiang, Zhenxing Ren, Ziwei Song, Yahui Zhang, Xingxing Guo, Genquan Han, and Yue Hao.  Computing primitive of fully vcsel-based all-optical spiking neural network for supervised learning and pattern classification.  IEEE Transactions on Neural Networks and Learning Systems, 32(6):2494–2505, 2020.
  * [12] Jinlong Xiang, Axel Torchy, Xuhan Guo, and Yikai Su.  All-optical spiking neuron based on passive microresonator.  Journal of Lightwave Technology, 38(15):4019–4029, 2020.
  * [13] Johannes Feldmann, Nathan Youngblood, C David Wright, Harish Bhaskaran, and Wolfram HP Pernice.  All-optical spiking neurosynaptic networks with self-learning capabilities.  Nature, 569(7755):208–214, 2019.
  * [14] Timothée Masquelier and Simon J Thorpe.  Unsupervised learning of visual features through spike timing dependent plasticity.  PLoS computational biology, 3(2):e31, 2007.
  * [15] Saeed Reza Kheradpisheh, Mohammad Ganjtabesh, Simon J Thorpe, and Timothée Masquelier.  Stdp-based spiking deep convolutional neural networks for object recognition.  Neural Networks, 99:56–67, 2018.
  * [16] Murat Yessenov, Basanta Bhaduri, Peter J Delfyett, and Ayman F Abouraddy.  Free-space optical delay line using space-time wave packets.  Nature Communications, 11(1):5782, 2020.
  * [17] Neil Zeghidour, Olivier Teboul, Félix de Chaumont Quitry, and Marco Tagliasacchi.  Leaf: A learnable frontend for audio classification.  arXiv preprint arXiv:2101.08596, 2021.
  * [18] Amirreza Ahmadnejad, Ahmad Mahmmodian Darviishani, Mohmmad Mehrdad Asadi, Sajjad Saffariyeh, Pedram Yousef, and Emad Fatemizadeh.  Tacnet: Temporal audio source counting network, 2023.
  * [19] Yoshua Bengio, Thomas Mesnard, Asja Fischer, Saizheng Zhang, and Yuhuai Wu.  Stdp-compatible approximation of backpropagation in an energy-based model.  Neural computation, 29(3):555–577, 2017.
  * [20] Judson P Jones and Larry A Palmer.  An evaluation of the two-dimensional gabor filter model of simple receptive fields in cat striate cortex.  Journal of neurophysiology, 58(6):1233–1258, 1987.
  * [21] Wolfgang Fuhl and Enkelejda Kasneci.  Multi layer neural networks as replacement for pooling operations.  arXiv preprint arXiv:2006.06969, 2020.
  * [22] Alex Krizhevsky, Ilya Sutskever, and Geoffrey E Hinton.  Imagenet classification with deep convolutional neural networks.  Communications of the ACM, 60(6):84–90, 2017.
  * [23] Bastian Leibe and Bernt Schiele.  Analyzing appearance and contour based methods for object categorization.  In 2003 IEEE Computer Society Conference on Computer Vision and Pattern Recognition, 2003. Proceedings., volume 2, pages II–409. IEEE, 2003.
  * [24] Li Fei-Fei, Robert Fergus, and Pietro Perona.  One-shot learning of object categories.  IEEE transactions on pattern analysis and machine intelligence, 28(4):594–611, 2006.
  * [25] Xianxin Guo, Thomas D Barrett, Zhiming M Wang, and AI Lvovsky.  Backpropagation through nonlinear units for the all-optical training of neural networks.  Photonics Research, 9(3):B71–B80, 2021.
  * [26] Xing Lin, Yair Rivenson, Nezih T Yardimci, Muhammed Veli, Yi Luo, Mona Jarrahi, and Aydogan Ozcan.  All-optical machine learning using diffractive deep neural networks.  Science, 361(6406):1004–1008, 2018.
  * [27] Shane Colburn, Yi Chu, Eli Shilzerman, and Arka Majumdar.  Optical frontend for a convolutional neural network.  Applied optics, 58(12):3179–3186, 2019.
  * [28] William B Levy and Victoria G Calvert.  Computation in the human cerebral cortex uses less than 0.2 watts yet this great expense is optimal when considering communication costs.  BioRxiv, pages 2020–04, 2020.
  * [29] George A Rakuljic and Victor Leyva.  Volume holographic narrow-band optical filter.  Optics letters, 18(6):459–461, 1993.
  * [30] Miguel Iglesias Olmedo, Tianjian Zuo, Jesper Bevensee Jensen, Qiwen Zhong, Xiaogeng Xu, Sergei Popov, and Idelfonso Tafur Monroy.  Multiband carrierless amplitude phase modulation for high capacity optical data links.  Journal of Lightwave Technology, 32(4):798–804, 2013.
  * [31] Alejandro Baldominos, Yago Saez, and Pedro Isasi.  A survey of handwritten character recognition with mnist and emnist.  Applied Sciences, 9(15):3169, 2019.
  * [32] OpenAI.  Chatgpt.  https://openai.com/research/chatgpt, September 2023.  Accessed: Date.


Polyelectrolytes, such as polyacrylic acid (PAA), can effectively mitigate CaCO3 scale formation.
Despite their success as antiscalants, the underlying mechanism of Ca^2+ binding to polyelectrolyte chains remains unresolved. Through all-atom molecular dynamics simulations, we construct an adsorption isotherm of Ca^2+ binding to sodium polyacrylate (NaPAA) and investigate the associated binding mechanism. We find that the number of calcium ions adsorbed [Ca^2+_ads] to the polymer saturates at moderately high concentrations of free calcium ions [Ca^2+_aq] in the solution. This saturation value is intricately connected with the binding modes accessible to Ca^2+ ions when they bind to the polyelectrolyte chain. We identify two dominant binding modes:
the first involves binding to at most two carboxylate oxygens on a polyacrylate chain, and the second, termed the high binding mode, involves binding to four or more carboxylate oxygens. As the concentration of free calcium ions [Ca^2+_aq] increases from low to moderate levels, the polyelectrolyte chain undergoes a conformational transition from an extended coil to a hairpin-like structure, enhancing the accessibility to the high binding mode. At moderate concentrations of [Ca^2+_aq], the high binding mode accounts for at least a third of all binding events. The chain's conformational change and its consequent access to the high binding mode is found to increase the overall Ca^2+ ion binding capacity of the polyelectrolyte chain.

§ INTRODUCTION

Divalent metal ions, such as Ca^2+, exhibit a pronounced affinity for dissolved anions like carbonate (CO3^2-) in an aqueous solution.

The ion pairs, which readily associate, precipitate out of solution and form solid deposits (scale) due to their low solubility limit. Scale formation presents challenges to residential and industrial piping systems, restricting the fluid flow and fouling components.[Antony et al., 2011, Tijing et al., 2010, Muryanto et al., 2014] These metal ions also form complexes with household products, such as detergents, and disrupt their efficacy.

Polyelectrolytes, such as polyacrylic acid (PAA), are commonly employed to mitigate scale formation.[Huang et al., 2007, Sinn et al., 2004, Yu et al., 2004, Aschauer et al., 2010, Backer et al., 2021] Although it is not fully understood what makes an effective antiscalant polyelectrolyte, a few mechanisms for scale prevention have been proposed. The polyelectrolytes could prevent nucleation via chelating metal ions from the solution.[Gindele et al., 2022, Sinn et al., 2004] The chelation reduces the concentration of metal ion and the likelihood of their association with the dissolved anions. Simultaneously, polyelectrolytes could also adsorb onto the surfaces of scale crystals and prevent further growth or deposition.[Reddy and Hoch, 2001, Yu et al., 2004, Aschauer et al., 2010, Sparks et al., 2015]

The solution behavior of polyelectrolytes in divalent salt solutions poses important challenges in designing polyelectrolytes with enhanced antiscalant activity, as the polyelectrolyte–ion complex itself can precipitate and lead to further scale deposition. [Schweins and Huber, 2001, Schweins et al., 2003, Michaeli, 1960, Muthukumar, 2017, Schweins et al., 2006, Yethiraj, 2008, Park et al., 2012, Mantha and Yethiraj, 2015, Chialvo and Simonson, 2005, Dobrynin, 2006, Prabhu, 2005, Duan and Wang, 2022] Cloud point measurements revealed instances of phase separation into a polymer-poor (supernatant) liquid and a polymer-rich liquid with addition of divalent ions.[Sabbagh and Delsanti, 2000] Boisvert et al. performed osmotic pressure measurements and showed that the solution behavior of polyelectrolytes in divalent salt solutions is primarily influenced by a proposed site-binding mechanism of divalent cations[Boisvert et al., 2002]. The site-binding mechanism facilitates bridging of non-adjacent repeat units by the divalent cations. Through a Fourier Transform Infrared (FTIR) dialysis technique, [Fantinel et al., 2004] Fantinel et al. identified monodentate, bidentate, and bridging modes when Ca^2+ ions bind to carboxylate groups of a polyacrylate chain. However, the relative importance of each of these binding modes—particularly that of the bridging mode—on the ability of polyelectrolyte to chelate Ca^2+ ion has not been elucidated. Furthermore, the influence of Ca^2+ ion concentration on these binding modes remains elusive.

Mean-field theories have shed light on the solution behavior of polyelectrolytes in divalent salt solutions.[Solis and de laCruz, 2001, Kundagrami and Muthukumar, 2008, Wittmer et al., 1995, de la Cruz et al., 1995, Ermoshkin and de la Cruz, 2003, Lee and Muthukumar, 2009, Deserno et al., 1999] In addition to establishing conditions for the precipitation behaviours, these theoretical models have predicted a large reduction in polymer size with addition of divalent ions,[Liu et al., 2003]  beyond what is expected from electrostatic screening. The chain collapse was primarily attributed to ion bridging between non-neighboring repeat units. Coarse-grained implicit solvent molecular simulations, employing generic models for polyelectrolyte chains and ions, have confirmed the chain collapse and attributed it to the bridging capability of divalent ions.[Zhou et al., 2016, Liu et al., 2003, Kłos and Pakula, 2005] All-atom molecular dynamics simulations have gone a step further by explicitly treating solvent molecules to investigate the molecular principles that govern the binding of divalent cations to the polyelectrolyte chain.[Molnar and Rieger, 2004, Yao et al., 2018, Sappidi and Natarajan, 2016, Patel et al., 2017, Tribello et al., 2009, Chung and Huang, 2012, Bulo et al., 2007] These investigations have each reported that the Ca^2+ ion is strongly coordinated with the polyelectrolyte chain, resulting in highly coiled conformations with a chain rigidity reminiscent of crystal-like structures. Due to the strong Ca^2+–polyelectrolyte interactions and the number of Ca^2+ ions binding to the polyelectrolyte chain, these models hint at the overcharging of the Ca^2+–polyelectrolyte complex. However, recent potentiometric titration assays[Gindele et al., 2022] suggest that only $\frac{1}{3}$ of the binding sites are occupied before a fully charged polyelectrolyte chain reaches its saturated value of Ca^2+-binding capacity.

The saturation value in the adsorption isotherm, which describes the maximum amount of Ca^2+ that can be bound to a polyelectrolyte chain, reflects the chelating capacity of the polyelectrolyte chain. However, the molecular principles that govern the corresponding adsorption isotherm have not been addressed. Specifically, the interplay between the site-binding nature of Ca^2+–polyelectrolyte interactions and the conformational transitions of the polyelectrolyte, aimed at enhancing both the chelating capacity and the solubility of the Ca^2+–polyelectrolyte complex, has not been explored.

In this study, we address the mechanism of Ca^2+ adsorption onto a polyacrylate chain. We construct an adsorption isotherm to describe the binding behaviors of Ca^2+ ions to polyacrylate. We then determine the different binding modes accessible for Ca^2+ binding to the chain and quantify their impact on the adsorption isotherm. In our follow-up manuscript,[44] we address questions related to Ca^2+ ion-mediated association between polyelectrolyte chains. The rest of the manuscript is organized as follows. In <Ref>, we present a Hamiltonian Replica Exchange Molecular Dynamics (HREMD) protocol to selectively bias Ca^2+-polyelectrolyte interactions and efficiently sample the configurational space. We then introduce a free energy perturbation approach coupled with molecular dynamics to compute the adsorption isotherm describing Ca^2+ binding to the polyelectrolyte chain. We discuss the results of these calculations in <Ref> and present the conclusions in <Ref>.

§ MODELS AND METHODS

We investigated the mechanism of calcium ion binding to polyacrylate chain using all-atom molecular dynamics (MD) simulations. Our simulation system consists of a single polyacrylic acid (PAA) chain with 32 repeat units, solvated in a cubic water box with an edge length of 12 nm.  All repeat units on the polymer are charged, consistent with the solution conditions for antiscalant activity(i.e., solution pH $\thicksim 10$). Sodium ions were added for electroneutrality. The average end-to-end distance of such a polymer was $\thicksim 5$ nm. We chose our box dimensions so that the polymer does not interact with its periodic image. In <Ref>, we report compositions of different systems studied in this work.

$\mathbf{Label}$     $\mathbf{Ca^{2+}}$     $\mathbf{Cl^{-}}$     $\mathbf{Na^{+}}$     $\mathbf{Water}$
0CaCl_2     0     0     32     56,448
4CaCl_2     4     8     32     56,436
8CaCl_2     8     16     32     56,424
16CaCl_2     16     32     32     56,400
32CaCl_2     32     64     32     56,352
64CaCl_2     64     128     32     56,256
96CaCl_2     96     192     32     56,160
128CaCl_2     128     256     32     56,064

Composition of the cubic simulation box with an edge length of 12 nm, containing a single sodium polyacrylate chain with 32 repeat units, and at various numbers of CaCl2 in water.

Force field choice and the importance of solvent electronic polarization:
We employed the Generalized AMBER Force Field (GAFF) with the SPC/E water model, which had been rigorously tested by Mintis et al. for modeling the polyacrylate chain. [Mintis and Mavrantzas, 2019] During our modeling of calcium ions, we observed that the “full” electrostatic charge force field parameters overestimated calcium ion binding to carboxylate groups on the polyacrylate chain, resulting in a charge inversion of the polymer chain inconsistent with experimental reports (see Supporting Information). Duboué-Dijon et al., who investigated calcium ion binding to insulin, reported that such an overestimation was due to inadequate treatment of electronic dielectric screening when using full charges on ions with non-polarizable forcefields.[Duboué-Dijon et al., 2018] The correct energetics of ion-pair formation could be captured by molecular models with explicit polarization.[Bedrov et al., 2019] However, a general-purpose force field with explicit polarization, tested to reproduce the properties of polyacrylates, was not readily available, and polarizeable force fields introduce large computational expense.

Alternatively, in our models, we included electronic dielectric screening by uniformly scaling the charges of all solute atoms. This approach, known as the electronic continuum charge correction (ECC), is a mean-field method that attempts to mimic charge-carrying species within an electronic dielectric continuum.[Duboué-Dijon et al., 2020, Leontyev and Stuchebrukhov, 2011] While the ECC scheme is a physically meaningful concept, it is primarily used as an ad hoc solution to incorporate electronic polarization into an otherwise non-polarizable force field. ECC schemes tend to fail in systems with a discontinuity in the high-frequency dielectric constant.[Kirby and Jungwirth, 2019] Moreover, force field parameterizations implicitly account for electronic polarization effects to some extent. Introducing additional ECC correction may overly compensate for electronic polarization effects, resulting in a significant underestimation of cohesive energy density and leading to unphysical solution behavior.[Cui et al., 2019] Nevertheless, the ECC scheme yields meaningful observations in common electrolyte systems where the high-frequency dielectric constant remains uniform throughout.

Within the scope of our study, the interactions of interest involve calcium ions and carboxylate groups on the polyacrylate chain. Biomolecular systems with the same carboxylate-calcium ion pair have shown success with ECC schemes. [Timr et al., 2018, Martinek et al., 2018, Kohagen et al., 2014, Kohagen et al., 2014] We particularly employ the ECC scheme reported by Jungwirth and coworkers, who provided parameters for modeling calcium and other ions in an aqueous solution.[Martinek et al., 2018] In their work, the authors uniformly scaled the charges on ions by a factor of 0.75. Lennard-Jones parameters describing dispersion interactions of corresponding ions were optimized to reproduce ab initio molecular dynamics results for ion-pairing and neutron scattering experiments. Such a parameter set accurately described the properties of calcium ions in an aqueous solution and their association with carboxylate groups on amino acids. We used the same parameter set to model calcium and other ions in our simulation system.

When modeling the electronic polarization effects of the polyacrylate chain in an aqueous solution, we scaled the partial charges of all the atoms on the polyacrylate chain by 0.75 and used Lennard-Jones parameters reported by Mintis et al. [Mintis and Mavrantzas, 2019] to describe dispersion interactions. Although our approach is a commonly accepted practice,[Kohagen et al., 2014] we emphasize that it is not rigorous. However, different properties of polymers strongly depend on their chain length, and finding the right strategy to optimize their dispersion interaction parameters is not obvious.

Since the “full-charge” parameters by Mintis et al. [Mintis and Mavrantzas, 2019] reproduce the structural properties of the polyacrylate chain in a salt-free solution, we validated the predictions of the “scaled-charge” model against the former. Even though we did not re-optimize the Lennard-Jones interaction parameters of the polyelectrolyte chain, the conformational flexibility of a polyelectrolyte chain in a salt-free solution did not change. Additionally, we observed identical distributions for structural properties when compared to the full-charge” model (see Supporting Information).

Simulation methodology: We used GROMACS 2022.3 [Abraham et al., 2015, Spoel et al., 2005, Berendsen et al., 1995] patched with PLUMED 2.8.1 [Bonomi et al., 2009, Tribello et al., 2014, PlumedConsortium, 2019] and employed the following protocol to conduct our molecular dynamics simulations of the systems reported in <Ref>. First, we minimized the energy of the system using the steepest descent algorithm until the maximum force on any atom in the system was smaller than 100 kJ mol$^{-1}$ nm$^{-1}$. Next, we equilibrated the system at constant NPT conditions with pressure and temperature set to 1 atm and 300 K, respectively. While fluctuations in the energy and box size minimized within a few nanoseconds of the equilibration run, it took tens of nanoseconds to equilibrate the chain structure. We used the Berendsen barostat[Berendsen et al., 1984] with the velocity-rescaling stochastic thermostat during the first 10 ns of the equilibration. For the remainder of the equilibration run, we switched to the Parrinello–Rahman barostat[Parrinello and Rahman, 1981] with the Nosé–Hoover chain thermostat for the accurate reproduction of the thermodynamic properties of the system.[Nosé, 1984, Hoover, 1985] Following the equilibration run, we conducted production MD simulations of these systems in the NVT ensemble using the Nosé–Hoover chain thermostat.

We employed the leap-frog time integration algorithm with a finite time step of 2 fs to integrate the equations of motion. Additionally, we utilized the LINCS constraint algorithm to convert all bonds with hydrogen atoms into constraints.[Hess et al., 1997] We applied periodic boundary conditions along all three spatial axes and used the particle mesh Ewald (PME) method with a minimum Fourier spacing of 0.12 nm to calculate the long-ranged electrostatic interactions.[Darden et al., 1993, Essmann et al., 1995] We applied a cutoff distance of 1.2 nm for computing van der Waals interactions. We used the same distance as the real-space cut-off value while computing the PME electrostatics.

Need for enhanced sampling molecular simulations: Although the ECC scheme with the non-polarizable force field greatly reduced the PAA–Ca^2+ binding/unbinding relaxation times, regular MD simulations were unable to efficiently sample the polymer conformational space in an aqueous CaCl2 solution. Even a microsecond long trajectory was not sufficient to sample polymer conformations in any of the systems listed in <Ref> with calcium numbers higher than 4 Ca ions. We direct readers to the Supporting Information document for the relevant data and discussion.

To address the challenges associated with polymer conformational sampling in an aqueous CaCl2 solution, we employed Hamiltonian Replica Exchange Molecular Dynamics (HREMD).[Wang et al., 2011, Bussi, 2013] Our HREMD framework is based on the flexible implementation of the REST2 variant,[Wang et al., 2011] as previously reported by Bussi.[Bussi, 2013] We introduced a parameter $\lambda$ to selectively bias the interactions between the polymer and ions, as well as the dihedral potential components of the Hamiltonian. The charges of the ions and the polymer were scaled by a factor of $\sqrt{\lambda}$, while their Lennard-Jones interaction parameter ($\epsilon$) was scaled by $\lambda$. Similarly, the polymer dihedral potential was also scaled by $\lambda$. With this scheme for $\lambda$-parameterized Hamiltonian, $\lambda=1$ corresponds to the system of interest with full-scale interactions. We determined that the parameterized Hamiltonian with $\lambda=0.67$ rapidly sampled the polymer conformational space.  Coordinate exchange between neighboring replicas were attempted every 500 steps. We utilized 16 replicas, with $\lambda$ values ranging from 1 to 0.67 (geometrically spaced), to simulate polymer conformations in an aqueous solution containing 8 CaCl2 or 16 CaCl2. For higher Ca$^{2+}$ numbers, we increased the number of replicas to 24. This combination of number of replicas and the $\lambda$-range yielded acceptable exchange probabilities ($\thicksim0.3$) between neighboring replicas. All the relevant results reported in the subsequent sections were obtained by averaging over a 250 ns production HREMD run, conducted at constant volume and a temperature of 300 K.

Computing ion adsorption isotherm from molecular dynamics simulations: The primary objective of the current work is to investigate Ca$^{2+}$ chelation onto a model polyacrylate chain. In an aqueous CaCl2 solution containing a polyacrylate chain, a dynamic equilibrium exists between calcium ions that are freely dispersed in the solution ($\mathrm{Ca^{2+}_{aq. free}}$) and the calcium ions adsorbed per monomer of the polyacrylate chain ($\mathrm{AA-Ca^{2+}}$).

\begin{equation}
\mathrm{AA} + \mathrm{Ca^{2+}_{aq.free}}\rightleftharpoons \mathrm{AA-Ca^{2+}}
\label{eq:Equil}
\end{equation}

Here, $\mathrm{AA}$ represents the concentration of repeat units on the polymer chain that are not bound to any calcium ions. An adsorption isotherm, quantifying <Ref>, describes calcium ion chelating ability of a model polyacrylate chain. We constructed the isotherm by plotting the number of calcium ions adsorbed per monomer ($\mathrm{AA-Ca^{2+}}$) against the concentration of calcium ions that are freely dispersed in the solution ($\mathrm{Ca^{2+}_{aq. free}}$).

Computing $\mathrm{AA-Ca^{2+}}$ from the simulation trajectory is straightforward. We calculated the number of calcium ions that are within 0.7 nm (see SI) of the polymer atoms at each frame. This quantity was then divided by the number of repeat units per chain (i.e. 32 in this case) and ensemble averaged over the trajectory.

We determined $\mathrm{Ca^{2+}_{aq. free}}$ by equating the chemical potential of the Ca$^{2+}$ ions in the system ($\mu^\mathrm{System}_{\ce{CaCl2}}$) with that of a pure CaCl$_2$ aqueous suspension ($\mu^\mathrm{Solution}_{\ce{CaCl2}}$).
This required conducting simulations of two separate sets of systems: one containing a polyacrylate chain to determine $\mu^\mathrm{System}_{\ce{CaCl2}}$, and the other without a polyacrylate chain for $\mu^\mathrm{Solution}_{\ce{CaCl2}}$.
We employed the procedure laid out by Panagiotopoulos and coworkers to determine the chemical potentials via free energy perturbation approach.[Panagiotopoulos, 2020, Saravi and Panagiotopoulos, 2022, Young and Panagiotopoulos, 2018, Young et al., 2019]

In brief, $\mu_{\ce{CaCl2}}$ represents the change in free energy when adding (or removing) a unit of CaCl2,  This is expressed in <Ref> as the sum of the corresponding ideal gas component ($\mu^{\mathrm{id}}_{\mathrm{CaCl_2}}$) and the residual component ($\mu^{\mathrm{R}}_{\mathrm{CaCl_2}}$).

\begin{align}
&\mu_{\ce{CaCl2}}=\mu^{\mathrm{id}}_{\ce{CaCl2}} + \mu^{\mathrm{R}}_{\ce{CaCl2}} \nonumber \\
&\mu^{\mathrm{id}}_{\ce{CaCl2}} = \mu^0_{\ce{Ca^{2+}}} + 2\mu^0_{\ce{Cl^{-}}} + 3RT\ln{\left(\frac{\sqrt[3]{4}N_{\ce{CaCl2}}k_BT}{{P^0\left<V\right>}}\right)} \nonumber \\
&\mu^{\mathrm{R}}_{\ce{CaCl2}} = \mu^\mathrm{R}_\mathrm{vdw} + \mu^\mathrm{R}_\mathrm{coul}
\label{eq:mu}
\end{align}

Here, $\mu^0_{\ce{Cl^{-}}}$ and $\mu^0_{\ce{Ca^{2+}}}$ are the respective chemical potentials of Cl^- and Ca^2+ at the reference pressure of 1 bar. We used the tabulated value of $\mu^0_{\ce{Cl^{-}}}$ from the NIST-JANAF thermochemical tables[M. W. Chase, 1998]. Moučka et al. estimated the value for $\mu^0_{\ce{Ca^{2+}}}$[Moučka et al., 2018], and we employed the same value in our calculations. $\mathrm{N_{CaCl_2}}$ represents the number of CaCl2 units in the simulation box, $k_B$ denotes the Boltzmann constant, $P^0=1,\mathrm{bar}$ is the reference pressure, and $\left<V\right>$ stands for the average volume obtained from NPT simulations of a system with a specific composition.

To calculate $\mu^{\mathrm{R}}_{\mathrm{CaCl_2}}$, we gradually decoupled  a Ca^2+ and two Cl^- from the system. The initial configuration for this study came from the most probable polymer conformation identified from the previously described enhanced sampling simulations. We then tagged a Ca^2+ ion that was adsorbed onto the polymer backbone. In a solution without the polyacrylate chain, a randomly selected Ca^2+ ion served as the tagged ion, and a comparison point between the two systems. Regardless of the system, two randomly chosen Cl^- ions were tagged.

The electrostatic interactions of the tagged Ca^2+ ion with other particles in the system were gradually turned off in 11 stages ($\phi=1,0.9,\ldots,0$). Here, $\phi=1$ represents the system with fully activated electrostatic interactions of the tagged calcium ion, while $\phi=0$ corresponds to complete deactivation. Each stage consisted of a series of molecular dynamics simulation steps, beginning with energy minimization, followed by 5 ns of equilibration and a 10 ns production simulation run at  a pressure of 1 bar and a temperature of 300 K. The output from the production step of the $i$th stage ($\phi_i$) served as the initial configuration for the MD simulation steps of the $(i+1)$th stage ($\phi_{i+1}$). The same methodology was then applied to deactivate the electrostatic interactions of each tagged chloride ion and also the van der Waals interactions of the tagged calcium ion and the two tagged chloride ions. We then used the Bennett Acceptance Ratio (BAR) [Bennett, 1976, Mester and Panagiotopoulos, 2015] approach, natively implemented in the GROMACS MD package[Abraham et al., 2015], to estimate the residual chemical potential ($\mu^{\mathrm{R}}_{\ce{CaCl2}}$) from these different stages of MD simulations.  Notably, given the position-dependent nature of the unbinding free energy of a Ca^2+ ion from the polymer backbone, we calculated the unbinding free energy for each adsorbed Ca^2+ ion individually and used their average to estimate $\mu^{\mathrm{R}}_{\ce{CaCl2}}$.

§ RESULTS AND DISCUSSION

The ability of a polyacrylate chain to sequester Ca^2+ ions is intricately tied to the bulk solution concentration of Ca^2+, denoted as $\left[\mathrm{Ca^{2+}_{{aq.free}}}\right]$. In <Ref>, we present the adsorption isotherm and describe the extent of Ca^2+ binding onto a 32-mer polyacrylate chain as a function of $\left[\mathrm{Ca^{2+}_{{aq.free}}}\right]$. At low Ca^2+ concentrations, most of the binding sites on the polyacrylate chain are available, and an increase in $\left[\mathrm{Ca^{2+}_{{aq.free}}}\right]$ results in a rapid increase in the number of Ca^2+ ions sequestered by the polyacrylate chain. However, at moderate concentrations, the accessible binding sites become occupied, and the polyacrylate chain saturates with Ca^2+.

Isotherm describing Ca^2+ ion binding to a polyacrylate chain with 32 repeat units. Error bars represent one standard deviation around the sample mean.  Note: The apparent positive slope of the last two data points is a consequence of large uncertainties in determining corresponding solution concentration of free Ca^2+ in the system. The differences in their vertical axis values are very minimal.

From the plateau in the adsorption isotherm, we note that the maximum binding capacity of a model polyacrylate chain is $0.40 \pm 0.05$ Ca^2+ ions per monomer, which aligns with recent potentiometric titration experiments[Gindele et al., 2022]. Although the adsorption isotherm bears some resemblance to a Langmuir model, we observe that the Ca^2+ binding does not obey the model assumptions. Notably, the conformational flexibility of the polyelectrolyte chain results in distinct chemical environments around each binding site, rendering the binding sites non-equivalent.

We demonstrate this phenomenon in a polyelectrolyte solution corresponding to $\mathrm{Ca^{2+}_{aq.free}=0.026 mol/kg}$. First, we identify the system configuration that corresponds to the polymer chain with the most probable radius of gyration (<Ref>). Then, we independently unbind each of the bound Ca^2+ ions. We employ the free energy perturbation approach described in <Ref> to compute the unbinding free energy, and report the corresponding values in <Ref>.

Label     $\mathrm{\Delta F \left[kJ/mol\right]}$

1     $814\pm 1$

2     $815\pm 1$

3     $814\pm 1$

4     $818\pm 2$

5     $818\pm 2$

6     $821\pm 3$

7     $821\pm 3$

8     $815\pm 1$

9     $828\pm 2$

10     $823\pm 1$

11     $817\pm 1$

1*Average     1*$818\pm 5$

(sfig:BindSt) Conformation of a 32-mer polyacrylate chain with the most probable radius of gyration in a solution corresponding to $\mathrm{\left[Ca^{2+}\right]_{aq.free}=0.026 ~ mol/kg}$
and (sfig:UnEnr)
free energies of unbinding a Ca^2+ ion that is adsorbed to the polyacrylate chain.

The unbinding free-energies of the 11 Ca^2+ ions broadly fall into two classes.
We categorize these binding sites into a “low” binding mode (${\le 818}$ kJ/mol) and “high” binding mode (${> 818}$ kJ/mol). From the binding sites depicted in <Ref> along with the accompanying free energies in <Ref>, we identify that the high binding mode is approximately 5–10 kJ/mol energetically more favorable than the low binding mode and facilitates a Ca^2+ ion-bridge between non-neighboring carboxylate groups on the polyacrylate chain.

The various binding modes and their coupling to polyelectrolyte conformation could impact the number of Ca^2+ ions sequestered. We determine a binding mode by tracking the number of carboxylate oxygen atoms that are around a Ca^2+ ion. From the radial distribution of carboxylate oxygen atoms around a Ca^2+ ion, we identify that their most probable separation is about 0.35 nm (see Supporting Information). We compute the probability of finding a varying number of carboxylate oxygen atoms within 0.35 nm from a Ca^2+ ion and report this as a function of $\mathrm{Ca^{2+}_{{aq.free}}}$ in <Ref>.

Probability of finding certain number of carboxylate oxygen atoms ($\mathrm{O_{Carboxylate}}$) on a 32-mer polyacrylate chain around a Ca^2+ ion.

We observe that the low binding mode, in which two carboxylate oxygen atoms bind to a given Ca^2+ ion, remains dominant at all concentrations of $\mathrm{Ca^{2+}_{{aq.free}}}$ studied. However, with an increase in $\mathrm{Ca^{2+}_{{aq.free}}}$, the high binding mode corresponding to ion-bridging becomes more favorable. Interestingly, when $\mathrm{Ca^{2+}_{{aq.free}}}$ exceeds $2.18 \times 10^{-2}$ mol/kg, this trend reverses: further increases in $\mathrm{Ca^{2+}_{{aq.free}}}$ promote the low binding mode once more. We hypothesize that this non-monotonic trend in the population of different binding modes arises from the enhanced electrostatic screening at higher Ca^2+ concentrations.

Nevertheless, at moderate and high concentrations of $[\mathrm{Ca^{2+}_{{aq.free}}}$, nearly $\tfrac{1}{3}$ of the binding events are due to the high binding mode. This high coordination binding environment is seen to have a large impact on the polymer size and conformation. We investigate the coupling between the binding modes and the polyelectrolyte chain conformation by tracking the chain radius of gyration ($\mathrm{R_g}$) as a function of $\mathrm{Ca^{2+}_{{aq.free}}}$. We report these results in <Ref>.

a) Radius of gyration of 32-mer polyacrylate chain at different concentrations of $\left[\mathrm{Ca^{2+}_{{aq.free}}}\right]$. Dominant binding modes at b) low, c) moderate, and d) high $\mathrm{\left[Ca^{2+}\right]_{aq.free}}$ concentrations. The conformation in (c) corresponds to the $\mathrm{R}_\mathrm{g}$ minimum in (a)

We note from <Ref>(a) that at low concentrations of $\mathrm{Ca^{2+}_{{aq.free}}}$, where the population of the high binding mode is insignificant, the polymer chain adopts an extended conformation with an $\mathrm{R_g}$ of approximately 1.65–1.9 nm (<Ref>(b)). At these concentrations, Ca^2+ ions bind to at most one carboxylate group on the polyacrylate chain. As $\mathrm{Ca^{2+}_{{aq.free}}}$ increases and the population of high coordination binding sites subsequently increases, the polymer conformation transitions to a collapsed state (<Ref>(c)). Here, the high coordination binding sites, which bridge two strands of the polyelectrolyte chain, are nearly as prominent as the low coordination binding sites located on the solvent-exposed side of each strand. Intriguingly, in solutions with high concentrations of Ca^2+ ions, although some of the bridging events are disrupted, the conformation of a polyelectrolyte chain still resembles that of a collapsed state (<Ref>(d)). The disruption of the bridging events due to enhanced screening from the other ions in the system increases the conformational flexibility and hence the average chain size.

Since the high coordination binding sites are energetically more favorable, we anticipate that access to a larger population of these binding sites would enhance the ability of a polyelectrolyte chain to sequester Ca^2+ ions. To explore this, we investigate the inverse problem; we limit the polyelectrolyte chain to only binding sites with low coordination environments and construct the adsorption isotherm. We achieve this by first identifying the polymer conformation that corresponds to the most probable $\mathrm{R}_\mathrm{g}$ in a system with no added Ca^2+ ions. Utilizing this chain conformation, with harmonic restraints imposed, we prepare polyelectrolyte solutions with added Ca^2+ ions. We compute the Ca^2+ adsorption isotherm from these systems and compare it to that obtained from unrestrained simulations reported earlier. We show these results in <Ref>.

Ca^2+ ion binding isotherm to a polyacrylate chain with 32 repeat units restrained to an extended coil conformation and unrestricted.

We find that at low concentrations, the number of Ca^2+ ions adsorbed on the polyacrylate chain in both the restrained and unrestrained simulations is indistinguishable. This is because, at lower concentrations of Ca^2+ ions in the solution, unrestrained polyelectrolyte chains can only access the low coordination binding sites. However, we see a noticeable difference at moderate and high concentrations, where the unrestrained polyacrylate chain favors the high coordination binding sites. The restrained simulations, which do not have access to high binding modes, consistently adsorb fewer Ca^2+ ions per chain when compared to the unrestrained simulations.

These observations indicate that a polyelectrolyte chain's ability to efficiently chelate Ca^2+ ions is impacted by its access to a large number of high coordination binding sites. These sites are energetically more favorable for ion binding, and as such, they provide a more stable environment for the ions, leading to a more effective sequestration. This insight could prove useful in applications where selective and efficient ion capture is paramount, such as in water treatment processes or biomedical applications.

§ SUMMARY

Using all-atom molecular simulations coupled with a free energy perturbation approach, we constructed an adsorption isotherm to describe the binding of Ca^2+ ions to a model polyacrylate chain. Analysis of the adsorption isotherm revealed that the per-monomer Ca^2+ ion binding capacity of a fully charged polyelectrolyte chain saturates at a value of $0.40\pm0.05$. This saturation value correlates with the binding modes accessible to Ca^2+ ions as they bind to the polyelectrolyte chain. Two predominant binding modes were identified: one mode involves Ca^2+ ions binding to at most two carboxylate oxygen atoms on a polyacrylate chain, and the other involves Ca^2+ ions binding to four or more carboxylate oxygen atoms.

The population of low binding mode sites remains high across all concentrations of Ca^2+ in the solution. Nevertheless, at least one-third of the binding events at moderate and high concentrations of Ca^2+ ions in the solution are defined by the high binding mode events. These binding events, while responsible for enhancing the Ca^2+ ion binding capacity of a polyelectrolyte chain, also lead to the collapse of the polyelectrolyte chain's conformation. The solution concentration of Ca^2+ ions corresponding to the conformational transition falls within the range close to the supernatant (polymer-poor) side of the phase behavior reported by Sabbagh et al.[Sabbagh and Delsanti, 2000] The adsorption isotherm constructed in this study suggests that, at these concentrations, the collapsed polyelectrolyte chain had attained the saturation value of the Ca^2+ binding capacity.

The collapsed polyelectrolyte chain, with all available binding sites saturated, may indicate the onset of putative phase separation. This observation carries important implications for the chelating capacity of a polyelectrolyte towards Ca^2+ ions, and consequently, for its antiscalant activity.  Yet, polyelectrolyte-divalent ion complex falling out of equilibrium is inherently a multi-chain problem. In our follow up work,[44] we tackle this problem and  establish the conditions under which a polyelectrolyte in divalent salt solution falls out of equilibrium.

This work was supported by the Dow Chemical Company through a University Partnership Initiative (UPI) grant. We benefited greatly from the discussions with our Dow collaborators: Thomas Kalantar, Christopher Tucker, Larisa Reyes, and Meng Jing. S.M. thanks Dr. Chris Balzer for discussions and useful comments.

In the supporting document, we first demonstrate overcharging of the polyacrylate-Ca^2+ complex when using non-polarizable force fields with full charges on ions. Then, we validate scaled charge force fields for modeling polyacrylate conformations. Following that, we present evidence for meaningful sampling of polyelectrolyte-Ca^2+ complexes in an aqueous solution using a non-polarizable force field with electronic continuum correction and Hamiltonian Replica-Exchange Molecular simulations. Finally, we discuss how we established chemical potential equivalence conditions and determined the concentration of free Ca^2+ ions in the solution that is in equilibrium with the system containing the polyacrylate chain.

[Antony et al., 2011]
Antony, A.; Low, J. H.; Gray, S.; Childress, A. E.; Le-Clech, P.; Leslie, G.
Scale formation and control in high pressure membrane water treatment
systems: A review. Journal of Membrane Science 2011,
383, 1–16, DOI: 10.1016/j.memsci.2011.08.054
[Tijing et al., 2010]
Tijing, L. D.; Kim, H. Y.; Lee, D. H.; Kim, C. S.; Cho, Y. I. Physical water
treatment using RF electric fields for the mitigation of CaCO3 fouling in
cooling water. International Journal of Heat and Mass Transfer
2010, 53, 1426–1437, DOI:
[Muryanto et al., 2014]
Muryanto, S.; Bayuseno, A.; Ma'mun, H.; Usamah, M.; Jotho, Calcium Carbonate
Scale Formation in Pipes: Effect of Flow Rates, Temperature, and Malic Acid
as Additives on the Mass and Morphology of the Scale. Procedia
Chemistry 2014, 9, 69–76, DOI:
[Huang et al., 2007]
Huang, S.-C.; Naka, K.; Chujo, Y. Effect of Molecular Weights of Poly(acrylic
acid) on Crystallization of Calcium Carbonate by the Delayed Addition Method.
Polymer Journal 2007, 40, 154–162, DOI:
[Sinn et al., 2004]
Sinn, C. G.; Dimova, R.; Antonietti, M. Isothermal Titration Calorimetry of the
Polyelectrolyte/Water Interaction and Binding of
Ca^2+: Effects Determining the Quality of Polymeric
Scale Inhibitors. Macromolecules 2004, 37, 3444–3450,
DOI: 10.1021/ma030550s
[Yu et al., 2004]
Yu, J.; Lei, M.; Cheng, B.; Zhao, X. Effects of PAA additive and temperature
on morphology of calcium carbonate particles. Journal of Solid State
Chemistry 2004, 177, 681–689, DOI:
[Aschauer et al., 2010]
Aschauer, U.; Spagnoli, D.; Bowen, P.; Parker, S. C. Growth modification of
seeded calcite using carboxylic acids: Atomistic simulations. Journal
of Colloid and Interface Science 2010, 346, 226–231, DOI:
[Backer et al., 2021]
Backer, S.; Creamer, M.; Pulukkody, R.; Ferrieux, S. Title Not Provided. U.S.
Patent Application, 2021
[Gindele et al., 2022]
Gindele, M. B.; Malaszuk, K. K.; Peter, C.; Gebauer, D. On the Binding
Mechanisms of Calcium Ions to Polycarboxylates: Effects of Molecular Weight,
Side Chain, and Backbone Chemistry. Langmuir 2022, 38,
14409–14421, DOI: 10.1021/acs.langmuir.2c01662
[Reddy and Hoch, 2001]
Reddy, M. M.; Hoch, A. R. Calcite Crystal Growth Rate Inhibition by
Polycarboxylic Acids. Journal of Colloid and Interface Science
2001, 235, 365–370, DOI: 10.1006/jcis.2000.7378
[Sparks et al., 2015]
Sparks, D. J.; Romero-González, M. E.; El-Taboni, E.; Freeman, C. L.;
Hall, S. A.; Kakonyi, G.; Swanson, L.; Banwart, S. A.; Harding, J. H.
Adsorption of poly acrylic acid onto the surface of calcite: an experimental
and simulation study. Physical Chemistry Chemical Physics
2015, 17, 27357–27365, DOI: 10.1039/c5cp00945f
[Schweins and Huber, 2001]
Schweins, R.; Huber, K. Collapse of sodium polyacrylate chains in calcium salt
solutions. The European Physical Journal E 2001, 5,
117–126, DOI: 10.1007/s101890170093
[Schweins et al., 2003]
Schweins, R.; Lindner, P.; Huber, K. Calcium Induced Shrinking of NaPA
Chains: A SANS Investigation of Single Chain Behavior.
Macromolecules 2003, 36, 9564–9573, DOI:
[Michaeli, 1960]
Michaeli, I. Ion binding and the formation of insoluble polymethacrylic salts.
Journal of Polymer Science 1960, 48, 291–299, DOI:
[Muthukumar, 2017]
Muthukumar, M. 50th Anniversary Perspective: A Perspective on Polyelectrolyte
Solutions. Macromolecules 2017, 50, 9528–9560, DOI:
[Schweins et al., 2006]
Schweins, R.; Goerigk, G.; Huber, K. Shrinking of anionic polyacrylate coils
induced by Ca^2+, Sr^2+ and Ba^2+: A combined light
scattering and ASAXS study. The European Physical Journal E
2006, 21, 99–110, DOI:
[Yethiraj, 2008]
Yethiraj, A. Liquid State Theory of Polyelectrolyte Solutions. The
Journal of Physical Chemistry B 2008, 113, 1539–1551, DOI:
[Park et al., 2012]
Park, S.; Zhu, X.; Yethiraj, A. Atomistic Simulations of Dilute Polyelectrolyte
Solutions. The Journal of Physical Chemistry B 2012,
116, 4319–4327, DOI: 10.1021/jp208138t
[Mantha and Yethiraj, 2015]
Mantha, S.; Yethiraj, A. Conformational Properties of Sodium
Polystyrenesulfonate in Water: Insights from a Coarse-Grained Model with
Explicit Solvent. The Journal of Physical Chemistry B 2015,
119, 11010–11018, DOI: 10.1021/acs.jpcb.5b01700
[Chialvo and Simonson, 2005]
Chialvo, A. A.; Simonson, J. M. Solvation Behavior of Short-Chain Polystyrene
Sulfonate in Aqueous Electrolyte Solutions: A Molecular
Dynamics Study. The Journal of Physical Chemistry B 2005,
109, 23031–23042, DOI: 10.1021/jp053512e
[Dobrynin, 2006]
Dobrynin, A. V. Effect of Counterion Condensation on Rigidity of Semiflexible
Polyelectrolytes. Macromolecules 2006, 39, 9519–9527,
DOI: 10.1021/ma061030a
[Prabhu, 2005]
Prabhu, V. M. Counterion structure and dynamics in polyelectrolyte solutions.
Current Opinion in Colloid & Interface Science 2005,
10, 2–8, DOI: 10.1016/j.cocis.2005.04.002
[Duan and Wang, 2022]
Duan, C.; Wang, R. Association of two polyelectrolytes in salt solutions.
Soft Matter 2022, 18, 6934–6941, DOI:
[Sabbagh and Delsanti, 2000]
Sabbagh, I.; Delsanti, M. Solubility of highly charged anionic polyelectrolytes
in presence of multivalent cations: Specific interaction effect. The
European Physical Journal E - Soft Matter 2000, 1, 75–86,
DOI: 10.1007/s101890050009
[Boisvert et al., 2002]
Boisvert, J.-P.; Malgat, A.; Pochard, I.; Daneault, C. Influence of the
counter-ion on the effective charge of polyacrylic acid in dilute condition.
Polymer 2002, 43, 141–148, DOI:
[Fantinel et al., 2004]
Fantinel, F.; Rieger, J.; Molnar, F.; Hübler, P. Complexation of
Polyacrylates by Ca^2+ Ions. Time-Resolved Studies Using Attenuated
Total Reflectance Fourier Transform Infrared Dialysis Spectroscopy.
Langmuir 2004, 20, 2539–2542, DOI:
[Solis and de laCruz, 2001]
Solis, F.; de laCruz, M. O. Flexible linear polyelectrolytes in multivalent
salt solutions: Solubility conditions. The European Physical Journal E
2001, 4, 143–152, DOI: 10.1007/s101890170123
[Kundagrami and Muthukumar, 2008]
Kundagrami, A.; Muthukumar, M. Theory of competitive counterion adsorption on
flexible polyelectrolytes: Divalent salts. The Journal of Chemical
Physics 2008, 128, DOI: 10.1063/1.2940199
[Wittmer et al., 1995]
Wittmer, J.; Johner, A.; Joanny, J. F. Precipitation of Polyelectrolytes in the
Presence of Multivalent Salts. Journal de Physique II 1995,
5, 635–654, DOI: 10.1051/jp2:1995154
[de la Cruz et al., 1995]
de la Cruz, M. O.; Belloni, L.; Delsanti, M.; Dalbiez, J. P.; Spalla, O.;
Drifford, M. Precipitation of highly charged polyelectrolyte solutions in the
presence of multivalent salts. The Journal of Chemical Physics
1995, 103, 5781–5791, DOI: 10.1063/1.470459
[Ermoshkin and de la Cruz, 2003]
Ermoshkin, A. V.; de la Cruz, M. O. Polyelectrolytes in the Presence of
Multivalent Ions: Gelation Versus Segregation. Physical Review Letters
2003, 90, DOI: 10.1103/physrevlett.90.125504
[Lee and Muthukumar, 2009]
Lee, C.-L.; Muthukumar, M. Phase behavior of polyelectrolyte solutions with
salt. The Journal of Chemical Physics 2009, 130, DOI:
[Deserno et al., 1999]
Deserno, M.; Holm, C.; May, S. Fraction of Condensed Counterions around a
Charged Rod: Comparison of Poisson-Boltzmann Theory and
Computer Simulations. Macromolecules 1999, 33,
199–206, DOI: 10.1021/ma990897o
[Liu et al., 2003]
Liu, S.; Ghosh, K.; Muthukumar, M. Polyelectrolyte solutions with added salt: A
simulation study. The Journal of Chemical Physics 2003,
119, 1813–1823, DOI: 10.1063/1.1580109
[Zhou et al., 2016]
Zhou, J.; Barz, M.; Schmid, F. Complex formation between polyelectrolytes and
oppositely charged oligoelectrolytes. The Journal of Chemical Physics
2016, 144, DOI: 10.1063/1.4947255
[Kłos and Pakula, 2005]
Kłos, J.; Pakula, T. Monte Carlo simulations of a polyelectrolyte chain with
added salt: Effect of temperature and salt valence. The Journal of
Chemical Physics 2005, 123, DOI:
[Molnar and Rieger, 2004]
Molnar, F.; Rieger, J. “Like-Charge
Attraction” between Anionic
Polyelectrolytes: Molecular Dynamics Simulations.
Langmuir 2004, 21, 786–789, DOI:
[Yao et al., 2018]
Yao, G.; Zhao, J.; Ramisetti, S. B.; Wen, D. Atomistic Molecular Dynamic
Simulation of Dilute Poly(acrylic acid) Solution: Effects of Simulation Size
Sensitivity and Ionic Strength. Industrial & Engineering Chemistry
Research 2018, 57, 17129–17141, DOI:
[Sappidi and Natarajan, 2016]
Sappidi, P.; Natarajan, U. Effect of salt valency and concentration on
structure and thermodynamic behavior of anionic polyelectrolyte
Na^+-polyethacrylate aqueous solution. Journal of Molecular
Modeling 2016, 22, DOI:
[Patel et al., 2017]
Patel, K. H.; Chockalingam, R.; Natarajan, U. Molecular dynamic simulations
study of the effect of salt valency on structure and thermodynamic solvation
behaviour of anionic polyacrylate PAA in aqueous solutions. Molecular
Simulation 2017, 43, 691–705, DOI:
[Tribello et al., 2009]
Tribello, G. A.; Liew, C.; Parrinello, M. Binding of Calcium and Carbonate to
Polyacrylates. The Journal of Physical Chemistry B 2009,
113, 7081–7085, DOI: 10.1021/jp900283d
[Chung and Huang, 2012]
Chung, Y.-T.; Huang, C.-I. Ion condensation behavior and dynamics of water
molecules surrounding the sodium poly(methacrylic acid) chain in water: A
molecular dynamics study. The Journal of Chemical Physics
2012, 136, DOI: 10.1063/1.3697477
[Bulo et al., 2007]
Bulo, R. E.; Donadio, D.; Laio, A.; Molnar, F.; Rieger, J.; Parrinello, M. Site
Binding of Ca^2+ Ions to Polyacrylates in Water: A Molecular Dynamics
Study of Coiling and Aggregation. Macromolecules 2007,
40, 3437–3442, DOI: 10.1021/ma062467l
[44]
Glisman, A.; Mantha, S.; Yu, D.; Wang, Z.-G. Multivalent ion mediated
polyelectrolyte association and structure. Manuscript under review

[Mintis and Mavrantzas, 2019]
Mintis, D. G.; Mavrantzas, V. G. Effect of pH and Molecular Length on the
Structure and Dynamics of Short Poly(acrylic acid) in Dilute Solution:
Detailed Molecular Dynamics Study. The Journal of Physical Chemistry B
2019, 123, 4204–4219, DOI:
[Duboué-Dijon et al., 2018]
Duboué-Dijon, E.; Delcroix, P.; Martinez-Seara, H.;
Hladílková, J.; Coufal, P.; Křížek, T.;
Jungwirth, P. Binding of Divalent Cations to Insulin: Capillary
Electrophoresis and Molecular Simulations. The Journal of Physical
Chemistry B 2018, 122, 5640–5648, DOI:
[Bedrov et al., 2019]
Bedrov, D.; Piquemal, J.-P.; Borodin, O.; MacKerell, A. D.; Roux, B.;
Schröder, C. Molecular Dynamics Simulations of Ionic Liquids and
Electrolytes Using Polarizable Force Fields. Chemical Reviews
2019, 119, 7940–7995, DOI:
[Duboué-Dijon et al., 2020]
Duboué-Dijon, E.; Javanainen, M.; Delcroix, P.; Jungwirth, P.;
Martinez-Seara, H. A practical guide to biologically relevant molecular
simulations with charge scaling for electronic polarization. The
Journal of Chemical Physics 2020, 153, DOI:
[Leontyev and Stuchebrukhov, 2011]
Leontyev, I.; Stuchebrukhov, A. Accounting for electronic polarization in
non-polarizable force fields. Physical Chemistry Chemical Physics
2011, 13, 2613, DOI: 10.1039/c0cp01971b
[Kirby and Jungwirth, 2019]
Kirby, B. J.; Jungwirth, P. Charge Scaling Manifesto: A Way of Reconciling the
Inherently Macroscopic and Microscopic Natures of Molecular Simulations.
The Journal of Physical Chemistry Letters 2019, 10,
7531–7536, DOI: 10.1021/acs.jpclett.9b02652
[Cui et al., 2019]
Cui, K.; Yethiraj, A.; Schmidt, J. R. Influence of Charge Scaling on the
Solvation Properties of Ionic Liquid Solutions. The Journal of Physical
Chemistry B 2019, 123, 9222–9229, DOI:
[Timr et al., 2018]
Timr, Š.; Kadlec, J.; Srb, P.; Ollila, O. H. S.; Jungwirth, P. Calcium
Sensing by Recoverin: Effect of Protein Conformation on Ion Affinity.
The Journal of Physical Chemistry Letters 2018, 9,
1613–1619, DOI: 10.1021/acs.jpclett.8b00495
[Martinek et al., 2018]
Martinek, T.; Duboué-Dijon, E.; Timr, Š.; Mason, P. E.;
Baxová, K.; Fischer, H. E.; Schmidt, B.; Pluhařová, E.;
Jungwirth, P. Calcium ions in aqueous solutions: Accurate force field
description aided by ab initio molecular dynamics and neutron scattering.
The Journal of Chemical Physics 2018, 148, DOI:
[Kohagen et al., 2014]
Kohagen, M.; Mason, P. E.; Jungwirth, P. Accurate Description of Calcium
Solvation in Concentrated Aqueous Solutions. The Journal of Physical
Chemistry B 2014, 118, 7902–7909, DOI:
[Kohagen et al., 2014]
Kohagen, M.; Lepšík, M.; Jungwirth, P. Calcium Binding to
Calmodulin by Molecular Dynamics with Effective Polarization. The
Journal of Physical Chemistry Letters 2014, 5, 3964–3969,
DOI: 10.1021/jz502099g
[Abraham et al., 2015]
Abraham, M. J.; Murtola, T.; Schulz, R.; Páll, S.; Smith, J. C.;
Hess, B.; Lindahl, E. GROMACS: High performance molecular simulations
through multi-level parallelism from laptops to supercomputers.
SoftwareX 2015, 1-2, 19–25, DOI:
[Spoel et al., 2005]
Spoel, D. V. D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.;
Berendsen, H. J. C. GROMACS: Fast, flexible, and free. Journal of
Computational Chemistry 2005, 26, 1701–1718, DOI:
[Berendsen et al., 1995]
Berendsen, H.; van der Spoel, D.; van Drunen, R. GROMACS: A message-passing
parallel molecular dynamics implementation. Computer Physics
Communications 1995, 91, 43–56, DOI:
[Bonomi et al., 2009]
Bonomi, M.; Branduardi, D.; Bussi, G.; Camilloni, C.; Provasi, D.; Raiteri, P.;
Donadio, D.; Marinelli, F.; Pietrucci, F.; Broglia, R. A.; Parrinello, M.
PLUMED: A portable plugin for free-energy calculations with molecular
dynamics. Computer Physics Communications 2009, 180,
1961–1972, DOI: 10.1016/j.cpc.2009.05.011
[Tribello et al., 2014]
Tribello, G. A.; Bonomi, M.; Branduardi, D.; Camilloni, C.; Bussi, G. PLUMED
2: New feathers for an old bird. Computer Physics Communications
2014, 185, 604–613, DOI:
[PlumedConsortium, 2019]
PlumedConsortium, Promoting transparency and reproducibility in enhanced
molecular simulations. Nature Methods 2019, 16,
670–673, DOI: 10.1038/s41592-019-0506-8
[Berendsen et al., 1984]
Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.;
Haak, J. R. Molecular dynamics with coupling to an external bath. The
Journal of Chemical Physics 1984, 81, 3684–3690, DOI:
[Parrinello and Rahman, 1981]
Parrinello, M.; Rahman, A. Polymorphic transitions in single crystals: A new
molecular dynamics method. Journal of Applied Physics 1981,
52, 7182–7190, DOI: 10.1063/1.328693
[Nosé, 1984]
Nosé, S. A unified formulation of the constant temperature molecular
dynamics methods. The Journal of Chemical Physics 1984,
81, 511–519, DOI: 10.1063/1.447334
[Hoover, 1985]
Hoover, W. G. Canonical dynamics: Equilibrium phase-space distributions.
Phys. Rev. A 1985, 31, 1695–1697, DOI:
[Hess et al., 1997]
Hess, B.; Bekker, H.; Berendsen, H. J. C.; Fraaije, J. G. E. M. LINCS: A
linear constraint solver for molecular simulations. Journal of
Computational Chemistry 1997, 18, 1463–1472
[Darden et al., 1993]
Darden, T.; York, D.; Pedersen, L. Particle mesh Ewald: An $N \log {N}$) method
for Ewald sums in large systems. The Journal of Chemical Physics
1993, 98, 10089–10092, DOI: 10.1063/1.464397
[Essmann et al., 1995]
Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G.
A smooth particle mesh Ewald method. The Journal of Chemical Physics
1995, 103, 8577–8593, DOI: 10.1063/1.470117
[Wang et al., 2011]
Wang, L.; Friesner, R. A.; Berne, B. J. Replica Exchange with Solute Scaling: A
More Efficient Version of Replica Exchange with Solute Tempering (REST2).
The Journal of Physical Chemistry B 2011, 115,
9431–9438, DOI: 10.1021/jp204407d
[Bussi, 2013]
Bussi, G. Hamiltonian replica exchange in GROMACS: a flexible implementation.
Molecular Physics 2013, 112, 379–384, DOI:
[Panagiotopoulos, 2020]
Panagiotopoulos, A. Z. Simulations of activities, solubilities, transport
properties, and nucleation rates for aqueous electrolyte solutions. The
Journal of Chemical Physics 2020, 153, DOI:
[Saravi and Panagiotopoulos, 2022]
Saravi, S. H.; Panagiotopoulos, A. Z. Activity Coefficients and Solubilities of
NaCl in Water–Methanol Solutions from Molecular Dynamics Simulations.
The Journal of Physical Chemistry B 2022, 126,
2891–2898, DOI: 10.1021/acs.jpcb.2c00813
[Young and Panagiotopoulos, 2018]
Young, J. M.; Panagiotopoulos, A. Z. System-Size Dependence of Electrolyte
Activity Coefficients in Molecular Simulations. The Journal of Physical
Chemistry B 2018, 122, 3330–3338, DOI:
[Young et al., 2019]
Young, J. M.; Tietz, C.; Panagiotopoulos, A. Z. Activity Coefficients and
Solubility of CaCl2 from Molecular Simulations. Journal of
Chemical & Engineering Data 2019, 65, 337–348, DOI:
[M. W. Chase, 1998]
M. W. Chase, NIST-JANAF Thermochemical Tables - 4 ed. 1998;
[Moučka et al., 2018]
Moučka, F.; Kolafa, J.; Lísal, M.; Smith, W. R. Chemical
potentials of alkaline earth metal halide aqueous electrolytes and solubility
of their hydrates by molecular simulation: Application to CaCl2,
antarcticite, and sinjarite. The Journal of Chemical Physics
2018, 148, DOI: 10.1063/1.5024212
[Bennett, 1976]
Bennett, C. H. Efficient estimation of free energy differences from Monte Carlo
data. Journal of Computational Physics 1976, 22,
245–268, DOI: 10.1016/0021-9991(76)90078-4
[Mester and Panagiotopoulos, 2015]
Mester, Z.; Panagiotopoulos, A. Z. Mean ionic activity coefficients in aqueous
NaCl solutions from molecular dynamics simulations. The Journal of
Chemical Physics 2015, 142, DOI: