# The role of receptor uniformity in multivalent binding

Xiuyang Xia  School of Chemistry, Chemical Engineering and Biotechnology,
Nanyang Technological University, 62 Nanyang Drive, Singapore 637459  Division
of Physics and Applied Physics, School of Physical and Mathematical
Sciences,Nanyang Technological University, 21 Nanyang Link, Singapore 637371
Ge Zhang<EMAIL_ADDRESS>Department of Physics, City University of Hong
Kong, Hong Kong, China  Massimo Pica Ciamarra  Division of Physics and Applied
Physics, School of Physical and Mathematical Sciences,Nanyang Technological
University, 21 Nanyang Link, Singapore 637371  Yang Jiao  Materials Science
and Engineering, Arizona State University, Tempe, AZ 85287, USA  Ran Ni
<EMAIL_ADDRESS>School of Chemistry, Chemical Engineering and Biotechnology,
Nanyang Technological University, 62 Nanyang Drive, Singapore 637459

###### Abstract

Multivalency is prevalent in various biological systems and applications due
to the superselectivity that arises from the cooperativity of multivalent
binding. Traditionally, it was thought that weaker individual binding would
improve the selectivity in multivalent targeting. Here using analytical mean
field theory and Monte Carlo simulations, we discover that for receptors that
are highly uniformly distributed, the highest selectivity occurs at an
intermediate binding energy and can be significantly greater than the weak
binding limit. This is caused by an exponential relationship between the bound
fraction and receptor concentration, which is influenced by both the strength
and combinatorial entropy of binding. Our findings not only provide new
guidelines for the rational design of biosensors using multivalent nano-
particles but also introduce a new perspective in understanding biological
processes involving multivalency.

multivalent nano-particle binding, superselectivity, hyperuniformity,
combinatorial entropy, Monte Carlo simulation

##  I Introduction

Multivalent interactions play a crucial role in a variety of biological
processes Mammen  _et al._ (1998); Boudreau and Jones (1999); Huskens (2006);
Fasting  _et al._ (2012); Satav  _et al._ (2015); Karimi  _et al._ (2018).
They provide an “on-off” binding at a threshold receptor density, creating a
biological barcode - targeting surfaces that have a receptor density above the
threshold while leaving others untouched. As a result, the multivalent binding
strategy is also widely used in many biorelated applications, particularly in
drug delivery Baker (2010); Bartlett  _et al._ (2007); Davis  _et al._ (2010);
Akhtar  _et al._ (2014); Koenig  _et al._ (2021) and bio-sensing Zeng  _et
al._ (2012); Zhou  _et al._ (2014); Cai (2014).

The Martinez-Veracoecha and Frenkel (MF) model provides a selectivity
parameter $\alpha={\rm d}\ln\theta/{\rm d}\ln n_{R}$ quantifying the
dependence of targeted adsorption $\theta$ on the receptor density $n_{R}$
Martinez-Veracoechea and Frenkel (2011). Generally, the maximum of the
selectivity parameter $\alpha_{\rm max}$, where the targeted adsorption grows
fastest, is defined as the selectivity employed to characterize the overall
selectivity indicating the onset of guest nanoparticle binding and clustering
Martinez-Veracoechea and Frenkel (2011); Albertazzi  _et al._ (2013);
Angioletti-Uberti (2017). If the selectivity $\alpha_{\max}>1$, the binding of
nano-particles is superselective, which is a signature of multivalent binding.
The MF model predicts that $\alpha_{\rm max}$ increases as the binding
strength becomes weaker when neglecting nonspecific interactions Martinez-
Veracoechea and Frenkel (2011); Curk  _et al._ (2017), which was also observed
in recent experimental systems including DNA coated colloids Scheepers  _et
al._ (2020); Linne  _et al._ (2021); Zhang  _et al._ (2018), multivalent
guest-host polymers Albertazzi  _et al._ (2013); Dubacheva  _et al._ (2014,
2015, 2019) and influenza virus particles Overeem  _et al._ (2020).

All studies mentioned above assume that the receptors grafted on the host
substrate follow the Poisson distribution considering they are spatially
uncorrelated. However, due to the complex environment on cell membranes, the
receptors are heterogeneously distributed and correlated Lingwood and Simons
(2010), of which the effect remains unknown. Additionally, recent
breakthroughs in DNA nanotechnol- ogy offer the possibility to precisely
design the spatial distribution of receptors on a substrate Bila  _et al._
(2022). Here we investigate how the uniformity of receptor distribution
affects the selectivity in multivalent nano-particle binding by focusing on
the hyperuniform, Poisson, and anti-hyperuniform distributions Torquato  _et
al._ (2021). We find that the more uniformly distributed receptors lead to
higher selectivity $\alpha_{\max}$, and intriguingly, the maximum selectivity
appears at certain intermediate binding energy for hyperuniform distributions,
which is qualitatively different from the Poisson distribution and anti-
hyperuniform distributions with $\alpha_{\max}$ approaching the upper bound at
the infinitely weak binding energy limit. Moreover, the highest selectivity
obtained for receptors of hyperuniform distributions can be significantly
larger than the upper bound in the Poisson and anti-hyperuniform
distributions, where the relatively large number fluctuation of receptors
masks the effect and causes the selectivity to increase monotonically with
decreasing the binding strength.

## Methodology

As shown in Fig. 1a, we consider that immobile receptors are grafted on a host
substrate. The nano-particles are controlled by an activity
$z=v_{0}\exp(\beta\mu)$ with $\mu$ the chemical potential of nano-particles
and $\beta=1/k_{B}T$, where $v_{0}$ is the volume that each particle can
explore when bound on the substrate, and $k_{B}$ and $T$ are the Boltzmann
constant and temperature of the system, respectively. Each nano-particle is
coated with $\kappa$ mobile ligands, which can bind to the receptors
reversibly with the binding free energy $f_{B}$. The binding free energy
$f_{B}$ is determined by both the equilibrium constant of ligand-receptor
binding in solvent $K_{a}$ and the configurational entropy penalty due to the
constraint of tethering $\Delta S_{\rm conf}$: $\beta f_{B}=-\log
K_{a}-k_{B}^{-1}\Delta S_{\rm{conf}}$ Leunissen  _et al._ (2010). Assuming
that the adsorption of each guest particle is independent, we divide the
substrate into $N_{\rm max}$ sites, each of which can bind with one guest
particle at most. The fraction of sites that are occupied by particles with at
least one bond formed is

$\theta(z,n_{R})=\frac{zq(n_{R})}{1+zq(n_{R})}.$ (1)

Using the unbound site as the reference state, the single-site bound state
partition function $q(n_{R})$ with $n_{R}$ receptors can be written as

$q\left(n_{R}\right)=\sum_{\lambda=1}^{\min\left(\kappa,n_{R}\right)}Q(\lambda,n_{R}),$
(2)

where $\lambda$ is the number of bonds formed and

$Q(\lambda,n_{R})=e^{-\lambda\beta
f_{B}}\frac{\kappa!n_{R}!}{(\kappa-\lambda)!\lambda!\left(n_{R}-\lambda\right)!}.$
(3)

Figure 1:  Multivalent nano-particle binding. (a) Schematic representation of
the prototypical multivalent adsorption model, in which the ligands (blue) on
the particles (white) can bind with the immobile receptors (pink) on the
substrate (grey) reversibly. (b) Illustration of the $\kappa$-$\mu VT$ Monte
Carlo simulation, in which some receptors (red) are bound with implicit
ligands on the nano-particles (blue) while the others (black) are unbound. In
the simulation, the bonds are implicit. (c) Part of typical snapshots of
receptors following various distributions. The global typical snapshots of
receptors can be found in SI.

Then the fraction of bound sites or adsorption is

$\langle\theta\rangle=\left\langle\frac{zq}{1+zq}\right\rangle_{\langle
n_{R}\rangle},$ (4)

where $\langle\cdot\rangle_{\langle n_{R}\rangle}$ calculates the average over
the receptor number distribution with the mathematical estimate $\langle
n_{R}\rangle$. The selectivity parameter is defined as

$\alpha=\frac{{\rm d}\ln\langle\theta\rangle}{{\rm d}\ln\langle
n_{R}\rangle}.$ (5)

Since the higher selectivity usually appears at small activity Martinez-
Veracoechea and Frenkel (2011); Wang and Dormidontova (2012); Phan  _et al._
(2023), the time for the substrate to exchange nano-particles with the
reservoir to reach equilibrium is very long, which makes the direct Monte
Carlo (MC) simulations in 3D systems expensive and inefficient. Here we
propose a $\kappa$-$\mu VT$ MC simulation method with implicit ligands and
bonds in 2D, which enables us to efficiently sample in the additional bond
number dimension (see Methods). As shown in Fig. 1b, we model the multivalent
nano-particles as hard disks of diameter $\sigma$ and volume $s_{\rm
hd}=\pi\sigma^{2}/4$ controlled by the chemical potential $\mu$. We assume
that one receptor can only bind with the particle covering it, i.e., the
center-to-center distance between the receptor and ligand is less than
$\sigma/2$. The total number of sites $N_{\rm max}=L^{2}/s_{\rm hd}$, and the
average number of receptors per site $\langle n_{R}\rangle=N_{R}s_{\rm
hd}/L^{2}$. The activity $z=s_{\rm hd}\exp(\beta\mu)/\Lambda^{2}$ with
$\Lambda$ the de Broglie wavelength. The distribution $p(n_{R})$ is
numerically sampled by the number of receptors within a 2D spherical window of
radius $\sigma/2$. One can see that the $\kappa$-$\mu VT$ MC simulation
essentially simulates a monolayer of nano-particles near the host substrate,
where nano-particles can bind with receptors on the substrate, and the system
exchanges nano-particles with a bulk (3D) reservoir of chemical potential
$\mu$ above. The advantage of the $\kappa$-$\mu VT$ model is that one does not
need to explicitly simulate the exchange of nano-particles between the host
substrate and bulk reservoir through diffusion, which could be very
computationally expensive at small activity.

## Results

###  I.1 Receptor uniformity enhances selectivity

We consider four different types of receptor distributions: anti-hyperuniform
(Anti-HU) distributions Torquato (2018), the Poisson distribution, stealthy
hyperuniform (SHU) distributions Florescu  _et al._ (2009); Torquato  _et al._
(2015) and a square lattice (Fig. S1). In equilibrium, Anti-HU can describe
systems close to a critical point, and SHU describes the disordered systems
with long range correlations Torquato (2018). All those distributions are
statistically homogeneous point processes and follow the central limit
theorem, i.e., they can be approximated by Gaussian distributions at the large
$\langle n_{R}\rangle$ limit Torquato  _et al._ (2021). The spatial uniformity
of a receptor distribution at given $\langle n_{R}\rangle$ can be
characterized by the relative local number variance
$\sigma^{2}_{n_{R}}/\langle n_{R}\rangle$. For the Poisson distribution in 2D,
$\sigma^{2}_{n_{R}}/\langle n_{R}\rangle=1$. For a perfect square lattice,
$\sigma^{2}_{n_{R}}/\langle n_{R}\rangle\sim\langle n_{R}\rangle^{-1/2}$,
which essentially implies that the square lattice is more uniform than the
Poisson distribution.

SHU distributions follow the same scaling with the square lattice. The
configurations are generated by minimizing
$\Phi(\mathbf{r}^{N})=\sum_{|\mathbf{k}|<K}S(\mathbf{k})$ using the limited-
memory BFGS algorithm Liu and Nocedal (1989), starting from a Poisson
configuration with number density $\rho=1$. Here $K=4\sqrt{\pi\chi}$ Torquato
_et al._ (2015) and $\chi=M(K)/[D(N_{R}-1)]$ denotes the relative fraction of
constrained degrees of freedom compared to the total degrees of freedom
$D(N_{R}-1)$ with $M(K)$ the number of independently constrained wave vectors
Torquato (2018). The prefactor of the distributions depends on the parameter
$\chi$ with the larger $\chi$ being more uniform or with smaller density
fluctuations.

On the contrary, in Anti-HU structures, $\sigma^{2}_{n_{R}}$ increases faster
than $\langle n_{R}\rangle$, and here we choose configurations that exhibit
$\sigma^{2}_{n_{R}}/\langle n_{R}\rangle\sim\langle n_{R}\rangle^{1/2}$. The
anti-hyperuniform configurations are generated using the algorithm detailed in
Ref. Torquato (2021). Specifically, we use the limited-memory BFGS algorithm
to minimize $\sum_{|\mathbf{k}|<K}[\langle S(\mathbf{k})\rangle-
S_{0}(\mathbf{k})]^{2}$, in which $K=10$, $\langle S(\mathbf{k})\rangle$ is
the average structure factor
$S(\mathbf{k})=|\sum_{j=1}^{N_{R}}\exp(-i\mathbf{k}\cdot\mathbf{r}_{j})|^{2}/N$
over $N_{c}=100$ configurations, and the targeted structure factor for various
$a$ is $S_{0}(\mathbf{k})=1+a\exp(-|\mathbf{k}|)/|\mathbf{k}|$. The prefactor
of the distributions depends on the parameter $a$ with the larger $a$ being
less uniform or with the larger density fluctuations. For each structure, we
individually generate 10 snapshots of $10^{4}$ receptors to sample the spatial
distribution and to be used in MC simulations.

Figure 2:  Receptor uniformity enhances selectivity. Average bound fraction
$\langle\theta\rangle$ (a) and selectivity parameter $\alpha$ (b) as a
function of $\langle n_{R}\rangle$ with $\kappa=4$, $\beta f_{B}=-2$ and
$\beta\mu=-10$ for various receptors distributions. The solid curves are the
theoretical predictions of Eqs. 4 and 5, and the symbols are obtained from
simulations.  Figure 3:  Achieving the highest selectivity by tuning binding
free energy. $\langle\theta\rangle$ and $\alpha$ as a function of $\langle
n_{R}\rangle$ for various binding free energy $\beta f_{B}$ for typical
receptors structures: Anti-HU $a=10$, Poisson, SHU $\chi=0.48$ and square
lattice. The solid curves are the theoretical predictions of Eqs. 4 and 5, and
the symbols are obtained from simulations. In all simulations, $\kappa=4$ and
$\beta\mu=-10$.

In Fig. 2, we plot the average bound fraction $\langle\theta\rangle$ and the
selectivity parameter $\alpha$ as functions of $\langle n_{R}\rangle$ for
various receptor distributions. One can see that Eqs. 4 and 5 agree
quantitatively with computer simulations when $\langle\theta\rangle<0.5$
(indicated by the dotted horizontal line), and at very large $\langle
n_{R}\rangle$, the theoretically predicted $\langle\theta\rangle$ is larger.
This discrepancy is due to the fact that the excluded volume effect between
nano-particles is not considered in the mean field theory, which overestimates
the adsorption at high density. When $\langle n_{R}\rangle$ is small, with
increasing $\langle n_{R}\rangle$, less uniform distributions lead to larger
value of $\langle\theta\rangle$. It is because that
$\langle\theta\rangle\approx z\langle q\rangle$ according to Eq. 4, where the
average bound state partition function over the receptor distribution $\langle
q\rangle$ is a convex function (see SI), hence $\langle q\rangle$ increases
with increasing the variance of the distribution $\sigma^{2}_{n_{R}}$. As
shown in Fig. 2, with increasing the uniformity, i.e., from Anti-HU to the
Poisson, SHU structures and square lattice, the selectivity $\alpha_{\max}$
increases monotonically, and for SHU structures and square lattice,
$\alpha_{\max}$ is even larger than $\kappa=4$. This is intriguing as it has
been accepted that weaker binding energy enhances selectivity, of which the
upper bound of $\alpha_{\max}$ is $\kappa$ Martinez-Veracoechea and Frenkel
(2011).

Figure 4:  Multivalent binding of nano-particles in the zero variance
scenario. (a) The probability for $\lambda$ bonds formed on bound guest nano-
particles $Q(\lambda)/q$ as a function of $n_{R}$ for (upper) $\beta f_{B}=4$
and (lower) $-4$. (b,c) The zero variance bound state partition function $q$
and zero variance selectivity at low activity limit ${\alpha_{\rm zv,0}}$ as a
function of $n_{R}$. Open symbols are the analytical results from Eq. 2.
Colored dotted lines are from Eq. 6. Colored dashed lines are (b)
$q=n_{R}e^{-\beta f_{B}}$ and (c) Eq. 7, respectively. Black dashed and dotted
lines indicate ${\alpha_{\rm zv,0}}=1$ and $\kappa$, respectively. (d)
$\theta/z$ (upper colored solid curves), $q$ (upper black curve) and
${\alpha_{\rm zv}}$ (lower colored curves), ${\alpha_{\rm zv,0}}$ (lower black
curve) as a function of $n_{R}$ with $\beta f_{B}=4$ for various activity $z$.
Symbols indicate the position of $n_{R}^{\rm s}$, which is the crosspoint of
$\theta/z$ and $1/z$ (upper colored dashed curves). In all calculations,
$\kappa=10$.

###  I.2 Achieve the highest selectivity by tuning binding energy

In Fig. 3, we plot $\langle\theta\rangle$ and $\alpha$ as functions of
$\langle n_{R}\rangle$ for binding energy from strong ($\beta f_{B}=-6$) to
weak ($\beta f_{B}=4$) of various receptor distributions. For the receptor
structures of the Poisson and Anti-HU distributions, $\alpha_{\rm max}$
increases monotonically with increasing $\beta f_{B}$, namely weaker binding
enhances selectivity, while for SHU structures with $\chi=0.48$ and square
lattice, $\alpha_{\rm max}$ reaches the maximum at about $\beta f_{B}=-4$. To
understand this, we start with the selectivity parameter ${\alpha_{\rm
zv}}=\mathrm{d}\ln\theta/\mathrm{d}\ln n_{r}$ in the zero variance scenario,
which is the uniform limit of receptors with $\sigma_{n_{R}}^{2}=0$. As the
superselective adsorption of nano-particles of interest mostly occurs at low
activity, according to Eq. 1, when $zq\to 0$, $\theta(n_{R})\approx zq(n_{R})$
and ${\alpha_{\rm zv}}\approx{\mathrm{d}\ln q}/{\mathrm{d}\ln n_{R}}$. We
define the selectivity parameter at the low activity limit as ${\alpha_{\rm
zv,0}}={\mathrm{d}\ln q}/{\mathrm{d}\ln n_{R}}$ 111Numerically, we calculate
the zero variance selectivity using the finite difference of $q$ or $\theta$
in the log-log scale: ${\alpha_{\rm
zv}}=[\log\theta(n_{R}+1)-\log\theta(n_{R}-1)]/[\log(n_{R}+1)-\log(n_{R}-1)]$
and ${\alpha_{\rm zv,0}}=[\log q(n_{R}+1)-\log
q(n_{R}-1)]/[\log(n_{R}+1)-\log(n_{R}-1)]$, respectively. We plot the
probability of forming $\lambda$ bonds on the guest nano-particle in the bound
state in Fig 4a. For weak binding $\beta f_{B}=4$ and small $n_{R}$, one can
see $q\left(n_{R}\right)\approx Q(\langle\lambda\rangle_{\rm bound},n_{R})$
with the most probable bond number $\langle\lambda\rangle_{\rm bound}\approx
1$ (yellow region in upper panel of Fig 4a). This implies that there is only
one bond formed for the particle bound on the substrate, and the ligands on
each particle cannot bind cooperatively. As shown in Fig. 4b, this leads to a
linear dependence $q\approx n_{R}\kappa e^{-\beta f_{B}}$ and ${\alpha_{\rm
zv,0}}\approx 1$ with no superselectivity. With increasing $n_{R}$,
$\langle\lambda\rangle_{\rm bound}$ approaches $\kappa$ at the large $n_{R}$
limit due to the restriction from the number of ligands available on nano-
particles, and this leads to a power-law dependence of $q$ on $n_{R}$ (see
SI):

$q\approx Q(\lambda=\kappa,n_{R})=e^{-\kappa\beta
f_{B}}\frac{n_{R}!}{(n_{R}-\kappa)!}\approx\left(n_{R}e^{-\beta
f_{B}}\right)^{\kappa},$ (6)

and ${\alpha_{\rm zv,0}}\approx\kappa$. Here, $\kappa$ ligands on each
particle bind with crowded receptors together, and the emergent combinatorial
entropy induces the power-law dependence.

Differently, at the strong binding $\beta f_{B}=-4$, although the power-law
dependence also appears at the large $n_{R}$ limit, in the small $n_{R}$
regime, i.e., $1<n_{R}<\kappa$, as shown in the lower panel in Fig. 4a, the
most probable bond number $\langle\lambda\rangle_{\rm bound}\approx n_{R}$. It
is because that the most probable bond number is limited by the number of
receptors on each site of the host substrate, rather than the number of
ligands on the nanoparticles, which implies an exponential dependence of $q$
on $n_{R}$ (see SI):

$q\approx Q(\lambda=n_{R},n_{R})=e^{-n_{R}\beta f_{B}}\frac{\kappa!}{(\kappa-
n_{R})!}\approx\left(\kappa e^{-\beta f_{B}}\right)^{n_{R}},$ (7)

and ${\alpha_{\rm zv,0}}\approx n_{R}(\ln\kappa-\beta f_{B})\sim n_{R}$. These
are the major results of this work. The derivation of Eqs. 6 and 7 based on
the saddle-point approximation method Angioletti-Uberti (2017); Varilly  _et
al._ (2012); Angioletti-Uberti  _et al._ (2013) can be found in SI. As shown
in Fig. 4c, ${\alpha_{\rm zv,0}}$ peaks around $n_{R}=\kappa$, which is
independent of the binding energy, and the selectivity ${\alpha_{\rm
zv,0}^{\rm max}}\approx\kappa(\ln\kappa-\beta f_{B})$ can be larger than
$\kappa$ when $\beta f_{B}<\ln\kappa-1$. Hence, in the zero variance binding
scenario, $\ln\kappa-1$ is a threshold of binding energy, lower than which
${\alpha_{\rm zv,0}^{\rm max}}>\kappa$, otherwise ${\alpha_{\rm zv,0}^{\rm
max}}$ has an upper bound $\kappa$ at $n_{R}\to+\infty$.

As shown in Fig. 4d, the rescaled bound fractions $\theta/z$ at various
activity $z$ collapse at small $n_{R}$, while reach reach the plateaus of
different height at high $n_{R}$. The behavior is qualitatively the same when
considering receptor distributions (Fig. S3). This indicates that $z$ does not
affect the selectivity parameter when $n_{R}<n_{R}^{\rm s}$, where
$q(n_{R}^{\rm s})\approx 1/z$ is the bound fraction saturation threshold.
Therefore, the zero variance selectivity parameter ${\alpha_{\rm
zv}}\approx{\alpha_{\rm zv,0}}$ at $n_{R}$ lower than the threshold, and drops
to 0 with further increasing $n_{R}$. Moreover, when $zq(n_{R}=1)>1$, i.e.,
$\beta f_{B}<\ln(z\kappa)$, the bound fraction saturates even if one nano-
particle only binds to one receptor, and no superselectivity occurs. To sum
up, an exponential dependence occurs at $\ln(z\kappa)<\beta f_{B}<\ln\kappa-1$
and $1<n_{R}<\kappa$ in the zero variance scenario.

For the highly uniformly distributed receptors, e.g., SHU $\chi=0.48$ and the
square lattice in Fig. 3, the situation is qualitatively the same to the zero
variance scenario since the variance of the receptor distribution is small.
The largest $\alpha_{\max}$ occurs at an intermediate binding energy with a
larger value than $\kappa$. Additionally, the critical $\langle n_{R}\rangle$
with $\alpha_{\max}$ appears independent of the binding energy $\beta
f_{B}=-6$ and $-4$, which is a distinct feature of the exponential dependence.

For the Poisson and less uniform receptor distributions, i.e., Anti-HU, our
numerical evidence shows that $\alpha_{\max}$ increases monotonically with
increasing $\beta f_{B}$ and approaches $\kappa$ at $\beta f_{B}\to+\infty$
(Fig. 3). We believe that this is due to the relatively large receptor number
fluctuations masking the exponential dependence in the intermediate binding
energy. As shown in Fig. S4, in the weak binding limit, receptor uniformity
has little effect on $\langle\theta\rangle$, and $\alpha_{\max}\to\kappa$
holds as long as the receptor distribution obeys the central limit theorem and
satisfies $\sigma_{n_{R}}^{2}<\langle n_{R}\rangle^{2}$ at large $\langle
n_{R}\rangle$ (see SI). Therefore, weaker binding enhances the selectivity
with an upper bound limit $\kappa$ for those less uniform distributions
Martinez-Veracoechea and Frenkel (2011).

Figure 5: Superselectivity of multivalent nano-particle binding on receptors
of tunable local uniformity. Selectivity $\alpha_{\max}$ as a function of
binding strength $\beta f_{B}$ for structures obtained from equilibrium hard-
disk fluids at various packing fraction $\phi$. In all simulations, $\kappa=4$
and $\beta\mu=-10$.

Next we investigate the binding of multivalent nano-particles to receptors
with tunable local uniformity. We use the configurations of receptors obtained
from equilibrium fluids of hard disks (HD) of various packing fraction $\phi$.
When $\phi\to 0$, the HD system recovers an ideal gas of the Poisson
distribution, i.e, $\sigma_{n_{R}}^{2}/\langle n_{R}\rangle=1$, and with
increasing $\phi$, at certain length scale, $\sigma_{n_{R}}^{2}/\langle
n_{R}\rangle\sim\langle n_{R}\rangle^{\xi}$ with $\xi<0$, which is the local
uniformity induced by the increased short range correlation because of the
excluded volume effect (Fig. S5). The local uniformity of the configuration
increases with increasing $\phi$ of HD systems, while at large enough length
scale $\sigma_{n_{R}}^{2}/\langle n_{R}\rangle\sim\langle n_{R}\rangle^{0}$
(Fig. S5). The measured $\alpha_{\max}$ as functions of binding strength
$\beta f_{B}$ for different $\phi$ are shown in Fig. 5. One can see that with
small local uniformity, e.g., $\phi=0.01$, $\alpha_{\max}$ increases
monotonically with increasing $\beta f_{B}$, which is the same as receptors of
the Poisson distribution (Fig. S6). However, with increasing $\phi$, the non-
monotonic dependence of $\alpha_{\max}$ on $\beta f_{B}$ appears when
$\phi>0.15$ (Fig. S6). This implies that the local uniformity induced by the
excluded volume effect can trigger the non-monotonic dependence of selectivity
on the binding free energy.

Here we consider the binding of multivalent nanoparticles on a rigid flat
substrate, while on a cell membrane, as more bonds form with the
nanoparticles, the membrane roughness decreases. Because the bonds suppress
membrane shape fluctuations, which promotes the formation of additional bonds
cooperatively Krobath  _et al._ (2009); Steinkühler  _et al._ (2019). In our
model, we can consider the thermal roughness of flexible substrates by
rewriting the partition function of $\lambda$ bonds forming in Eq. 3 to
include an entropy cost, $\Delta S_{\mathrm{mem}}(\lambda)$, that originates
from the suppression of membrane shape fluctuations upon $\lambda$ bonds
formation. This term depends on the relative roughness of the membranes,
$\xi_{\perp}\sim\lambda^{1/2}$, and a characteristic length $\xi_{RL}$ that
represents the extension of the receptor-ligand complex perpendicular to the
membranes Hu  _et al._ (2013). When $\xi_{\perp}\gg\xi_{RL}$, we have
$e^{k_{B}^{-1}\Delta
S_{\mathrm{mem}}(\lambda)}\sim\xi_{\perp}^{-1}\sim\lambda^{0.5}$ Krobath  _et
al._ (2007). At the large $n_{R}$ limit, where $q\approx
Q\left(\lambda=\kappa,n_{R}\right)$, the entropy cost term $\Delta
S_{\mathrm{mem}}$ is a constant because of the constant number of bonds, and
Eq. 6 remains valid. At this limit, the suppression of membrane shape
fluctuation has no qualitative effect on $\alpha_{zv,0}$. However, an
additional term should be added to $\alpha_{zv,0}$ to account for the
contribution of $\Delta S_{\mathrm{mem}}$. This term is given by ${{\rm d}\log
e^{k_{B}^{-1}\Delta S_{\mathrm{mem}}\left(n_{R}\right)}}/{{\rm d}\log
n_{R}}\approx 0.5$, which suggests that the suppression of membrane shape
fluctuations can enhance the superselectivity at this limit.

##  II Discussion

In conclusion, we have investigated the impact of receptor uniformity on the
superselective binding of multivalent nano-particles, for which we devised a
$\kappa$-$\mu VT$ MC simulation method to compare with analytical theory
without any fitting parameter. We find that receptors that are more uniformly
distributed lead to stronger superselective binding of multivalent particles.
Specifically, for receptors with SHU structures and square lattice, the
selectivity, $\alpha_{\max}$, can be significantly larger than the valence of
the nano-particle, $\kappa$, which is the highest level of selectivity that
receptors with Poisson and Anti-HU distributions can achieve. Furthermore, for
receptors with SHU distributions and square lattice arrangements, the largest
$\alpha_{\max}$ occurs at an intermediate strength of binding energy, which is
due to the exponential dependence of the bound fraction on the receptor
density. The exponential dependence arises from the restriction of available
receptors and is affected by both the binding energy and combinatorial
entropy. This is different from the binding on receptors with Poisson or Anti-
HU distributions, where weaker binding always enhances the selectivity. These
results suggest that for receptors that are highly uniformly distributed, one
does not have to use very weak binding energy to achieve high selectivity, and
the largest $\alpha_{\max}$ occurs at a certain relatively strong binding that
is less affected by non-specific bindings.

Our findings are relevant for designing superselective sensors using
multivalent nano-particles, of which a possible design is shown in Fig. 6.
Based on our results, one can use arrays of orderly arranged receptors of
different density using DNA origami to enhance the superselectivity using
strongly binding multivalent nano-particles to avoid the influence of non-
specific bindings Bila  _et al._ (2022). This also suggests the possibility of
designing superselective assembly of receptor-patterned colloidal systems
Mognetti  _et al._ (2019); Song  _et al._ (2020). Additionally, our findings
emphasize the significance of receptor distribution in biological systems.
Although in principle, many receptors on cell membranes are mobile, the
receptor diffusion time scale on cell membranes, $t_{\rm diff}$, could be much
longer than the time it takes to form/break a bond with ligands on nano-
particles, $t_{\rm on/off}$. This suggests that they can be considered as
effectively immobile receptors. Furthermore, there are also many immobile
receptors on cell membranes that are restricted and compartmentalized due to
interactions with the cytoskeleton Choquet and Triller (2003). Moreover, we
also show that the local uniformity induced by excluded volume effects can
trigger the non-monotonic dependence of selectivity on the binding free
energy, which suggests that our finding can be expected in systems of
relatively densely packed receptors, like the situation in many biological
membranes.

Figure 6: A possible design of superselective sensor for multivalent nano-
particles. The design consists of arrays of orderly arranged receptors of
density $n_{R}$ increasing from left to right, which can be realized by using
DNA origami Bila  _et al._ (2022). The inset shows typical adsorption curves
$\theta$ of multivalent nano-particles on the bands of receptors of different
$n_{R}$ grafted on the substrate.

##  III Supporting Information Available

  * •

SI.pdf: Simulation methods, additional theory details and figures of the main
text.

##  IV acknowledgments

We would like to thank Profs. Daan Frenkel and Stefano Angioletti-Uberti for
helpful discussions. This work is supported by the Academic Research Fund from
Singapore Ministry of Education Tier 1 Gant (RG59/21) and Tier 2 Grant
(MOE2019-T2-2-010).

## References

## References

  * Mammen  _et al._ (1998) M. Mammen, S.-K. Choi, and G.M. Whitesides, “Polyvalent interactions in biological systems: implications for design and use of multivalent ligands and inhibitors,” Angewandte Chemie International Edition 37, 2754–2794 (1998).
  * Boudreau and Jones (1999) N.J. Boudreau and P. L. Jones, “Extracellular matrix and integrin signalling: the shape of things to come,” Biochemical Journal 339, 481–488 (1999).
  * Huskens (2006) J. Huskens, “Multivalent interactions at interfaces,” Current opinion in chemical biology 10, 537–543 (2006).
  * Fasting  _et al._ (2012) C. Fasting, C.A. Schalley, M. Weber, O. Seitz, S. Hecht, B. Koksch, J. Dernedde, C. Graf, E.-W. Knapp, and R. Haag, “Multivalency as a chemical organization and action principle,” Angewandte Chemie International Edition 51, 10472–10498 (2012).
  * Satav  _et al._ (2015) T. Satav, J. Huskens, and P. Jonkheijm, “Effects of variations in ligand density on cell signaling,” Small 11, 5184–5199 (2015).
  * Karimi  _et al._ (2018) F. Karimi, A.J. O’Connor, G.G. Qiao, and D.E. Heath, “Integrin clustering matters: a review of biomaterials functionalized with multivalent integrin-binding ligands to improve cell adhesion, migration, differentiation, angiogenesis, and biomedical device integration,” Advanced healthcare materials 7, 1701324 (2018).
  * Baker (2010) M. Baker, “Homing in on delivery,” Nature 464, 1225–1227 (2010).
  * Bartlett  _et al._ (2007) D.W. Bartlett, H. Su, I.J. Hildebrandt, W.A. Weber, and M.E. Davis, “Impact of tumor-specific targeting on the biodistribution and efficacy of sirna nanoparticles measured by multimodality in vivo imaging,” Proceedings of the National Academy of Sciences 104, 15549–15554 (2007).
  * Davis  _et al._ (2010) M.E. Davis, J.E. Zuckerman, C.H.J Choi, D. Seligson, A. Tolcher, C.A. Alabi, Y. Yen, J.D. Heidel, and A. Ribas, “Evidence of rnai in humans from systemically administered sirna via targeted nanoparticles,” Nature 464, 1067–1070 (2010).
  * Akhtar  _et al._ (2014) M.J. Akhtar, M. Ahamed, H.A. Alhadlaq, S.A. Alrokayan, and S. Kumar, “Targeted anticancer therapy: overexpressed receptors and nanotechnology,” Clinica chimica acta 436, 78–92 (2014).
  * Koenig  _et al._ (2021) P.-A. Koenig, H. Das, H. Liu, B.M. Kümmerer, F.N. Gohr, L.-M. Jenster, L.D.J. Schiffelers, Y.M. Tesfamariam, M. Uchima, J.D. Wuerth, _et al._ , “Structure-guided multivalent nanobodies block sars-cov-2 infection and suppress mutational escape,” Science 371, eabe6230 (2021).
  * Zeng  _et al._ (2012) X. Zeng, C.A.S. Andrade, M.D.L. Oliveira, and X.-L. Sun, “Carbohydrate–protein interactions and their biosensing applications,” Analytical and bioanalytical chemistry 402, 3161–3176 (2012).
  * Zhou  _et al._ (2014) G. Zhou, M. Lin, P. Song, X. Chen, J. Chao, L. Wang, Q. Huang, W. Huang, C. Fan, and X. Zuo, “Multivalent capture and detection of cancer cells with dna nanostructured biosensors and multibranched hybridization chain reaction amplification,” Analytical chemistry 86, 7843–7848 (2014).
  * Cai (2014) W. Cai,  _Engineering in translational medicine_ (Springer, 2014).
  * Martinez-Veracoechea and Frenkel (2011) F.J. Martinez-Veracoechea and D. Frenkel, “Designing super selectivity in multivalent nano-particle binding,” Proceedings of the National Academy of Sciences 108, 10963–10968 (2011).
  * Albertazzi  _et al._ (2013) L. Albertazzi, F.J. Martinez-Veracoechea, C.M.A. Leenders, I.K. Voets, D. Frenkel, and E.W. Meijer, “Spatiotemporal control and superselectivity in supramolecular polymers using multivalency,” Proceedings of the National Academy of Sciences 110, 12203–12208 (2013).
  * Angioletti-Uberti (2017) S. Angioletti-Uberti, “Exploiting receptor competition to enhance nanoparticle binding selectivity,” Physical Review Letters 118, 068001 (2017).
  * Curk  _et al._ (2017) T. Curk, J. Dobnikar, and D. Frenkel, “Optimal multivalent targeting of membranes with many distinct receptors,” Proceedings of the National Academy of Sciences 114, 7210–7215 (2017).
  * Scheepers  _et al._ (2020) M.R.W. Scheepers, L.J. van IJzendoorn, and M.W.J. Prins, “Multivalent weak interactions enhance selectivity of interparticle binding,” Proceedings of the National Academy of Sciences 117, 22690–22697 (2020).
  * Linne  _et al._ (2021) C. Linne, D. Visco, S. Angioletti-Uberti, L. Laan, and D.J. Kraft, “Direct visualization of superselective colloid-surface binding mediated by multivalent interactions,” Proceedings of the National Academy of Sciences 118 (2021).
  * Zhang  _et al._ (2018) Y. Zhang, X. He, R. Zhuo, R. Sha, J. Brujic, N. Seeman, and P. Chaikin, “Multivalent, multiflavored droplets by design,” Proc. Natl Acad. Sci. USA 115, 9086–9091 (2018).
  * Dubacheva  _et al._ (2014) G.V. Dubacheva, T. Curk, B.M. Mognetti, R. Auzély-Velty, D. Frenkel, and R.P. Richter, “Superselective targeting using multivalent polymers,” Journal of the American Chemical Society 136, 1722–1725 (2014).
  * Dubacheva  _et al._ (2015) G.V. Dubacheva, T. Curk, R. Auzély-Velty, D. Frenkel, and R.P. Richter, “Designing multivalent probes for tunable superselective targeting,” Proceedings of the National Academy of Sciences 112, 5579–5584 (2015).
  * Dubacheva  _et al._ (2019) G.V. Dubacheva, T. Curk, D. Frenkel, and R.P. Richter, “Multivalent recognition at fluid surfaces: The interplay of receptor clustering and superselectivity,” Journal of the American Chemical Society 141, 2577–2588 (2019).
  * Overeem  _et al._ (2020) N.J. Overeem, P.H.E. Hamming, O.C. Grant, D. Di Iorio, M. Tieke, M.C. Bertolino, Z. Li, G. Vos, R.P. De Vries, R.J. Woods, _et al._ , “Hierarchical multivalent effects control influenza host specificity,” ACS central science 6, 2311–2318 (2020).
  * Lingwood and Simons (2010) D. Lingwood and K. Simons, “Lipid rafts as a membrane-organizing principle,” science 327, 46–50 (2010).
  * Bila  _et al._ (2022) H. Bila, K. Paloja, V. Caroprese, A. Kononenko, and M. Bastings, “Multivalent pattern recognition through control of nano-spacing in low-valency super-selective materials,” J. Am. Chem. Soc. 144, 21576 (2022).
  * Torquato  _et al._ (2021) S. Torquato, J. Kim, and M.A. Klatt, “Local number fluctuations in hyperuniform and nonhyperuniform systems: Higher-order moments and distribution functions,” Physical Review X 11, 021028 (2021).
  * Leunissen  _et al._ (2010) Mirjam E Leunissen, Remi Dreyfus, Roujie Sha, Nadrian C Seeman, and Paul M Chaikin, “Quantitative study of the association thermodynamics and kinetics of dna-coated particles for different functionalization schemes,” Journal of the American Chemical Society 132, 1903–1913 (2010).
  * Wang and Dormidontova (2012) S. Wang and E.E. Dormidontova, “Selectivity of ligand-receptor interactions between nanoparticle and cell surfaces,” Physical review letters 109, 238102 (2012).
  * Phan  _et al._ (2023) Huu Trong Phan, Dominic Lauzon, Alexis Vallée-Bélisle, Stefano Angioletti-Uberti, Jeanne Leblond Chain, and Suzanne Giasson, “Bimodal brush-functionalized nanoparticles selective to receptor surface density,” Proceedings of the National Academy of Sciences 120, e2208377120 (2023).
  * Torquato (2018) S. Torquato, “Hyperuniform states of matter,” Physics Reports 745, 1–95 (2018).
  * Florescu  _et al._ (2009) M. Florescu, S. Torquato, and P.J. Steinhardt, “Designer disordered materials with large, complete photonic band gaps,” Proceedings of the National Academy of Sciences 106, 20658–20663 (2009).
  * Torquato  _et al._ (2015) S. Torquato, G. Zhang, and F.H. Stillinger, “Ensemble theory for stealthy hyperuniform disordered ground states,” Physical Review X 5, 021020 (2015).
  * Liu and Nocedal (1989) Dong C Liu and Jorge Nocedal, “On the limited memory bfgs method for large scale optimization,” Mathematical programming 45, 503–528 (1989).
  * Torquato (2021) S. Torquato, “Structural characterization of many-particle systems on approach to hyperuniform states,” Physical Review E 103, 052126 (2021).
  * Note (1) Numerically, we calculate the zero variance selectivity using the finite difference of $q$ or $\theta$ in the log-log scale: ${\alpha_{\rm zv}}=[\mathop{log}\nolimits\theta(n_{R}+1)-\mathop{log}\nolimits\theta(n_{R}-1)]/[\mathop{log}\nolimits(n_{R}+1)-\mathop{log}\nolimits(n_{R}-1)]$ and ${\alpha_{\rm zv,0}}=[\mathop{log}\nolimits q(n_{R}+1)-\mathop{log}\nolimits q(n_{R}-1)]/[\mathop{log}\nolimits(n_{R}+1)-\mathop{log}\nolimits(n_{R}-1)]$, respectively.
  * Varilly  _et al._ (2012) P. Varilly, S. Angioletti-Uberti, B.M. Mognetti, and D. Frenkel, “A general theory of dna-mediated and other valence-limited colloidal interactions,” J. Chem. Phys. 137, 094108 (2012).
  * Angioletti-Uberti  _et al._ (2013) S. Angioletti-Uberti, P. Varilly, B.M. Mognetti, A.V. Tkachenko, and D. Frenkel, “Communication: A simple analytical formula for the free energy of ligand–receptor-mediated interactions,” J. Chem. Phys. 138, 01B401 (2013).
  * Krobath  _et al._ (2009) Heinrich Krobath, Bartosz Różycki, Reinhard Lipowsky, and Thomas R Weikl, “Binding cooperativity of membrane adhesion receptors,” Soft Matter 5, 3354–3361 (2009).
  * Steinkühler  _et al._ (2019) Jan Steinkühler, Bartosz Różycki, Cory Alvey, Reinhard Lipowsky, Thomas R Weikl, Rumiana Dimova, and Dennis E Discher, “Membrane fluctuations and acidosis regulate cooperative binding of ‘marker of self’protein cd47 with the macrophage checkpoint receptor sirp$\alpha$,” Journal of Cell Science 132, jcs216770 (2019).
  * Hu  _et al._ (2013) Jinglei Hu, Reinhard Lipowsky, and Thomas R Weikl, “Binding constants of membrane-anchored receptors and ligands depend strongly on the nanoscale roughness of membranes,” Proceedings of the National Academy of Sciences 110, 15283–15288 (2013).
  * Krobath  _et al._ (2007) HGJS Krobath, GJ Schütz, R Lipowsky, and TR Weikl, “Lateral diffusion of receptor-ligand bonds in membrane adhesion zones: Effect of thermal membrane roughness,” Europhysics Letters 78, 38003 (2007).
  * Mognetti  _et al._ (2019) Bortolo Matteo Mognetti, Pietro Cicuta, and Lorenzo Di Michele, “Programmable interactions with biomimetic dna linkers at fluid membranes and interfaces,” Reports on progress in physics 82, 116601 (2019).
  * Song  _et al._ (2020) Jake Song, Mehedi H Rizvi, Brian B Lynch, Jan Ilavsky, David Mankus, Joseph B Tracy, Gareth H McKinley, and Niels Holten-Andersen, “Programmable anisotropy and percolation in supramolecular patchy particle gels,” ACS nano 14, 17018–17027 (2020).
  * Choquet and Triller (2003) D. Choquet and A. Triller, “The role of receptor diffusion in the organization of the postsynaptic membrane,” Nature Reviews Neuroscience 4, 251–265 (2003).

Graphical TOC Entry


# Skin Lesion Classification Using Hybrid Deep Neural Networks

Amirreza Mahbod<EMAIL_ADDRESS>Institute of
Pathophysiology and Allergy Research, Medical University of Vienna, Vienna,
Austria  TissueGnostics GmbH, Vienna, Austria  Rupert Ecker  TissueGnostics
GmbH, Vienna, Austria  Isabella Ellinger  Institute of Pathophysiology and
Allergy Research, Medical University of Vienna, Vienna, Austria

###### Abstract

Skin cancer is one of the major types of cancers and its incidence has been
increasing over the past decades. Skin lesions can arise from various
dermatologic disorders and can be classified to various types according to
their texture, structure, color and other morphological features. The accuracy
of diagnosis of skin lesions, specifically the discrimination of benign and
malignant lesions, is paramount to ensure appropriate patient treatment.
Machine learning-based classification approaches are among popular automatic
methods for skin lesion classification. While there are many existing methods,
convolutional neural networks (CNN) have shown to be superior over other
classical machine learning methods for object detection and classification
tasks.
In this work, a fully automatic computerized method is proposed, which employs
well established pre-trained convolutional neural networks and ensembles
learning to classify skin lesions. We trained the networks using 2000 skin
lesion images available from the ISIC 2017 challenge, which has three main
categories and includes 374 melanoma, 254 seborrheic keratosis and 1372 benign
nevi images. The trained classifier was then tested on 150 unlabeled images.
The results, evaluated by the challenge organizer and based on the area under
the receiver operating characteristic curve (AUC), were 84.8% and 93.6% for
Melanoma and seborrheic keratosis binary classification problem, respectively.
The proposed method achieved competitive results to experienced dermatologist.
Further improvement and optimization of the proposed method with a larger
training dataset could lead to a more precise, reliable and robust method for
skin lesion classification.

##  1 Introduction

Skin cancer is the most common cancer type in the United States of America and
it is estimated that one in five Americans will develop skin cancer in their
lifetime. Among different types of skin cancer, malignant melanoma itself (the
deadliest type) is responsible for 10000 deaths annually in the united states
[1]. However, in case of diagnosis in early stages, it is completely curable
[2]. Early diagnosis of the cancer type is also important since different
treatments should be applied based on the type of cancer[3].
Dermoscopy has become the most important non-invasive tool for detection of
melanoma and other pigmented skinlesions. In contrast to naked eye
examinations, dermoscopy reveals subsurface structures of the skin. It makes
use of handheld devices that extend optical light ray penetration beyond the
skin surface and minimize surface reflection. [4]. This technique provides
better differentiation between different lesion types based on their
appearance and morphological features [5]. However, visual inspection of
dermoscopic images by dermatologists relies on the experience of the
physicians. Due to the difficulty and subjectivity of human interpretation,
computerized analysis of dermoscopy images has become an important research
area with the aim to reduce diagnostic errors [6]. While many computer-aided
methods for lesion classification are based on lesion border detection and
hand-crafted features, they lacked generalization capability due to high
variations in dermoscopic images. These variations could be caused by
different zooming, angles and lightening conditions [2],[7]. In addition,
dermoscopic images suffer from various artifacts, which could have negative
effect on classification accuracy. Fig. 1 shows some common artifacts.

Figure 1: Common artifacts in dermoscopic images. First row: a normal image,
dark corner artifact and skin hair artifact from left to right, respectively.
Second row: color chart artifact, ruler marker artifact and ink marker
artifact from left to right, respectively. (Adapted from ISIC archive)

These problems could be overcome by using established pre-trained deep
convolutional neural networks (CNN), which were trained on more that 1.4
million natural images (ImageNet dataset [8]) with various photographic
conditions. These networks can be used as optimized feature extractors for
dermoscopic images. While some research has been conducted to extract features
using deep features for skin lesion classification [3],[9], they were limited
by exploiting specific network architecture or using specific layer for
extracting features. In this work, we hypothesize that using different
architectures, extracting features from different layers and using ensemble
learning on top of all, could lead to better classification with competitive
precision to experienced dermatologist.

##  2 Method

The following subsections describe the algorithms and the dataset, which was
used in the frame of this paper.

###  2.1 Dataset

The available images from ISIC 2017 challenge were used for training purposes.
This dataset is composed of 2000 color dermoscopic skin images with
corresponding labels. Three different skin lesion types are found in this
dataset including 374 melanoma images, 254 seborrheic keratosis and 1372
benign nevi images with various sizes (from $1022\times 767$ to $6748\times
4499$), photographic angles and lightening conditions. Some images from this
dataset include artefacts as shown in Fig. 1. Another set of 150 skin images
was provided by the challenge organizers. The class labels of these test
images were concealed from the participants for final evaluation of the
performance of each proposed method by the challenge organizers.

###  2.2 Preprocessing

We applied three preprocessing steps on the images. First, we normalized the
images by subtracting the mean RGB value of ImageNet dataset as suggested in
[10] since the pre-trained networks were optimized for those images. Next, the
images were resized to appropriate size using bicubic interpolation to be fed
to the networks ($227\times 227$ and $224\times 224$). Finally, we augmented
the images to have a bigger training set. For data augmentation we rotated the
images by 0, 90, 180 and 270 degree and then applied horizontal flipping to
each image. From the derived artificial data set, we used 80% of the images
randomly for training the classifiers.

###  2.3 Deep learning models and classification

The presented deep feature extractor used two pre-trained networks. We used
pre-trained AlexNet architecture [10] and VGG-16 architecture [11] as feature
extractor. We explore extracting features from different layers of pre-trained
models to see how they could affect classification results. The features were
mainly extracted from the last fully connected (FC) layers. We used first and
second fully connected layers (referred as FC6 and FC7 with 4096 dimensions)
and also the concept detector layer (referred as FC8 with 1000 dimension).
These features along with corresponding labels (skin type) were then used to
train multi-class non-linear support vector machine (SVM) classifiers. We
trained different classifiers for each network and then averaged the class
score to have the final classification result.

###  2.4 Evaluation

To evaluate the classification results, we mapped SVM scores to probabilities
using logistic regression as described in [12]. Since the classifiers were
trained for a multi classification problem (3 classes including melanoma,
seborrheic keratosis and benign nevi classes), we combined the scores to have
the results for two binary classification problems, which were melanoma
classification and seborrheic keratosis classification. Classifier
performance, in terms of accuracy (ACC) and area under the receiver operating
characteristics curve (AUC) were measured for 150 test images.

##  3 Results

The summarized classification results on the test images are shown in Table 1.
To calculate the accuracy, a probability of 50% was used as threshold. All
experiments were performed on a laptop computer with NVIDIA GTX 840M 2GB , 6
GB installed memory and Intel Core i7-4500U 1.8 GHz CPU. The main
implementation of the algorithm was done with MATLAB (version 2016b) and
MatConvNet framework for CNN. The run-time of the implemented algorithm for
training SVM and extracting features was around 70 minutes and for testing on
all test images it was around 5 minutes.

Table 1:  Summarized results. In the table, M stands for Melanoma and SK
stands for seborrheic keratosis.

Network |  Feature Layer |  M-ACC |  SK-ACC |  M-AUC |  SK-AUC |  Avg AUC
---|---|---|---|---|---|---
AlexNet |  FC8 |  83.3 |  88.0 |  78.4 |  90.3 |  84.3
ALexNet |  All FC |  84 |  89.3 |  82.8 |  94.9 |  88.9
VGG-16 |  FC8 |  82.7 |  79.3 |  78.8 |  81 |  79.9
VGG-16 |  All FC |  85.3 |  88.0 |  81.9 |  92.8 |  87.3
Fusion |  All FC |  84.7 |  84.7 |  84.8 |  93.6 |  89.2

##  4 Discussion & Conclusion

In this work, we demonstrated how using pre-trained deep architectures, which
were originally trained for natural images, can be used for dermoscopic image
classification with competitive accuracy to experienced dermotologists. As the
results of Table 1 suggested, fusing networks results and extracting features
from different layers of the network could increase both accuracy and AUC.
Extracting features from other pre-trained networks and from other layers is
an open field for further research. In addition, using the same approach on a
larger training data set could potentially lead to a better predictive model.

## References

  * [1] H. W. Rogers, M. A. Weinstock, S. R. Feldman, and B. M. Coldiron, “Incidence estimate of nonmelanoma skin cancer (keratinocyte carcinomas) in the US population, 2012,” _JAMA dermatology_ , vol. 151, no. 10, pp. 1081–1086, 2015.
  * [2] G. Schaefer, B. Krawczyk, M. E. Celebi, and H. Iyatomi, “An ensemble classification approach for melanoma diagnosis,” _Memetic Computing_ , vol. 6, no. 4, pp. 233–240, 2014.
  * [3] J. Kawahara, A. BenTaieb, and G. Hamarneh, “Deep features to classify skin lesions,” in _Biomedical Imaging (ISBI), 2016 IEEE 13th International Symposium on_. IEEE, 2016, pp. 1397–1400.
  * [4] G. Argenziano, H. P. Soyer, V. De Giorgi, D. Piccolo, P. Carli, and M. Delfino, “Dermoscopy: a tutorial,” _EDRA, Medical Publishing & New Media_, p. 35, 2002.
  * [5] A. Steiner, M. Binder, M. Schemper, K. Wolff, and H. Pehamberger, “Statistical evaluation of epiluminescence microscopy criteria for melanocytic pigmented skin lesions,” _Journal of the American Academy of Dermatology_ , vol. 29, no. 4, pp. 581–588, 1993.
  * [6] M. G. Fleming, C. Steger, J. Zhang, J. Gao, A. B. Cognetta, and C. R. Dyer, “Techniques for a structural analysis of dermatoscopic imagery,” _Computerized medical imaging and graphics_ , vol. 22, no. 5, pp. 375–389, 1998.
  * [7] A. Esteva, B. Kuprel, R. A. Novoa, J. Ko, S. M. Swetter, H. M. Blau, and S. Thrun, “Dermatologist-level classification of skin cancer with deep neural networks,” _Nature_ , vol. 542, no. 7639, pp. 115–118, 2017.
  * [8] J. Deng, W. Dong, R. Socher, L.-J. Li, K. Li, and L. Fei-Fei, “Imagenet: A large-scale hierarchical image database,” in _Computer Vision and Pattern Recognition, 2009. CVPR 2009. IEEE Conference on_. IEEE, 2009, pp. 248–255.
  * [9] N. Codella, J. Cai, M. Abedini, R. Garnavi, A. Halpern, and J. R. Smith, “Deep learning, sparse coding, and SVM for melanoma recognition in dermoscopy images,” in _International Workshop on Machine Learning in Medical Imaging_. Springer, 2015, pp. 118–126.
  * [10] A. Krizhevsky, I. Sutskever, and G. E. Hinton, “Imagenet classification with deep convolutional neural networks,” in _Advances in neural information processing systems_ , 2012, pp. 1097–1105.
  * [11] K. Simonyan and A. Zisserman, “Very deep convolutional networks for large-scale image recognition,” _arXiv preprint arXiv:1409.1556_ , 2014\.
  * [12] J. Platt, “Probabilistic outputs for support vector machines and comparisons to regularized likelihood methods,” _Advances in large margin classifiers_ , vol. 10, no. 3, pp. 61–74, 1999.


# Magnetic Connectivity between Active Regions 10987, 10988, and 10989 by
Means of Nonlinear Force-Free Field Extrapolation

Tilaye Tadesse1 2 1 Max Planck Institut für Sonnensystemforschung, Max-Planck
Str. 2, D–37191 Katlenburg-Lindau, Germany email<EMAIL_ADDRESS>email<EMAIL_ADDRESS>email<EMAIL_ADDRESS>
2Addis Ababa University, College of Education, Department of Physics
Education, Po.Box 1176, Addis Ababa, Ethiopia
3 National Solar Observatory, Sunspot, NM 88349, U.S.A. email:
<EMAIL_ADDRESS>
T. Wiegelmann1 1 Max Planck Institut für Sonnensystemforschung, Max-Planck
Str. 2, D–37191 Katlenburg-Lindau, Germany email<EMAIL_ADDRESS>email<EMAIL_ADDRESS>email<EMAIL_ADDRESS>
2Addis Ababa University, College of Education, Department of Physics
Education, Po.Box 1176, Addis Ababa, Ethiopia
3 National Solar Observatory, Sunspot, NM 88349, U.S.A. email:
<EMAIL_ADDRESS>
and B. Inhester1 1 Max Planck Institut für Sonnensystemforschung, Max-Planck
Str. 2, D–37191 Katlenburg-Lindau, Germany email<EMAIL_ADDRESS>email<EMAIL_ADDRESS>email<EMAIL_ADDRESS>
2Addis Ababa University, College of Education, Department of Physics
Education, Po.Box 1176, Addis Ababa, Ethiopia
3 National Solar Observatory, Sunspot, NM 88349, U.S.A. email:
<EMAIL_ADDRESS>
and A. Pevtsov3 1 Max Planck Institut für Sonnensystemforschung, Max-Planck
Str. 2, D–37191 Katlenburg-Lindau, Germany email<EMAIL_ADDRESS>email<EMAIL_ADDRESS>email<EMAIL_ADDRESS>
2Addis Ababa University, College of Education, Department of Physics
Education, Po.Box 1176, Addis Ababa, Ethiopia
3 National Solar Observatory, Sunspot, NM 88349, U.S.A. email:
<EMAIL_ADDRESS>

###### Abstract

Extrapolation codes for modelling the magnetic field in the corona in
cartesian geometry do not take the curvature of the Sun’s surface into account
and can only be applied to relatively small areas, e.g., a single active
region. We apply a method for nonlinear force-free coronal magnetic field
modelling of photospheric vector magnetograms in spherical geometry which
allows us to study the connectivity between multi-active regions. We use
vector magnetograph data from the Synoptic Optical Long-term Investigations of
the Sun survey (SOLIS)/Vector Spectromagnetograph(VSM) to model the coronal
magnetic field, where we study three neighbouring magnetically connected
active regions (ARs: 10987, 10988, 10989) observed on 28, 29, and 30 March
2008, respectively. We compare the magnetic field topologies and the magnetic
energy densities and study the connectivities between the active regions(ARs).
We have studied the time evolution of magnetic field over the period of three
days and found no major changes in topologies as there was no major eruption
event. From this study we have concluded that active regions are much more
connected magnetically than the electric current.

###### keywords:

Magnetic fields Photosphere Corona

##  1 Introduction

S-Introduction In order to model and understand the physical mechanisms
underlying the various activity phenomena that can be observed in the solar
atmosphere, like, for instance, the onset of flares and coronal mass
ejections, the stability of active region, and to monitor the magnetic
helicity and free magnetic energy, the magnetic field vector throughout the
atmosphere must be known. However, routine measurements of the solar magnetic
field are mainly carried out in the photosphere. The magnetic field in the
photosphere is measured using the Zeeman effect of magnetically-sensitive
solar spectral lines. The problem of measuring the coronal field and its
embedded electrical currents thus leads us to use numerical modelling to infer
the field strength in the higher layers of the solar atmosphere from the
measured photospheric field. Except in eruptions, the magnetic field in the
solar corona evolves slowly as it responds to changes in the surface field,
implying that the electromagnetic Lorentz forces in this low-$\beta$
environment are relatively weak and that any electrical currents that exist
must be essentially parallel or antiparallel to the magnetic field wherever
the field is not negligible.

Due to the low value of the plasma $\beta$ (the ratio of gas pressure to
magnetic pressure), the solar corona is magnetically dominated [Gary (2001)].
To describe the equilibrium structure of the static coronal magnetic field
when non-magnetic forces are negligible, the force-free assumption is
appropriate:

$(\nabla\times\textbf{B})\times\textbf{B}=0\ilabel{one}$ (1)
$\nabla\cdot\textbf{B}=0\ilabel{two}$ (2)
$\textbf{B}=\textbf{B}_{\textrm{obs}}\quad\mbox{on photosphere}\ilabel{three}$
(3)

where B is the magnetic field and $\textbf{B}_{\textrm{obs}}$ is measured
vector field on the photosphere. Equation (one) states that the Lorentz force
vanishes (as a consequence of $\textbf{J}\parallel\textbf{B}$, where J is the
electric current density) and Equation (two) describes the absence of magnetic
monopoles.

The extrapolation methods based on this assumption are termed nonlinear force-
free field extrapolation [Sakurai (1981), Amari et al. (1997), Amari,
Boulmezaoud, and Mikic (1999), Amari, Boulmezaoud, and Aly (2006), Wu et al.
(1990), Cuperman, Demoulin, and Semel (1991), Demoulin, Cuperman, and Semel
(1992), Inhester and Wiegelmann (2006), Mikic and McClymont (1994),
Roumeliotis (1996), Yan and Sakurai (2000), Valori, Kliem, and Keppens (2005),
Wiegelmann (2004), Wheatland (2004), Wheatland and Régnier (2009), Wheatland
and Leka (2010), Amari and Aly (2010)]. For a more complete review of existing
methods for computing nonlinear force-free coronal magnetic fields, we refer
to the review papers by Amari97, Schrijver06, Metcalf, and Wiegelmann08.
Wiegelmann:2006 has developed a code for the self-consistent computation of
the coronal magnetic fields and the coronal plasma that uses non-force-free
MHD equilibria.

The magnetic field is not force-free in the photosphere, but becomes force-
free roughly 400 km above the photosphere [Metcalf et al. (1995)].
Furthermore, measurement errors, in particular for the transverse field
components (i.e., perpendicular to the line of sight of the observer), would
destroy the compatibility of a magnetogram with the condition of being force-
free. One way to ease these problems is to preprocess the magnetogram data as
suggested by Wiegelmann06sak. The preprocessing modifies the boundary values
of B within the error margins of the measurement in such a way that the moduli
of force-free integral constraints of Molodenskii [Molodensky (1974)] are
minimized. The resulting boundary values are expected to be more suitable for
an extrapolation into a force-free field than the original values.

In the present work, we use a larger computational domain which accommodates
most of the connectivity within the coronal region. We also take the
uncertainties of measurements in vector magnetograms into account as suggested
in DeRosa. We apply a preprocessing procedure to SOLIS data in spherical
geometry [Tadesse, Wiegelmann, and Inhester (2009)] by taking account of the
curvature of the Sun’s surface. For our observations, performed on 28, 29, and
30 March 2008, respectively, the large field of view contains three active
regions (ARs: 10987, 10988, 10989).

The full inversion of SOLIS/VSM magnetograms yields the magnetic filling
factor for each pixel, and it also corrects for magneto-optical effects in the
spectral line formation. The full inversion is performed in the framework of
Milne-Eddington model (ME)[Unno (1956)] only for pixels whose polarization is
above a selected threshold. Pixels with polarization below threshold are left
undetermined. These data gaps represent a major difficulty for existing
magnetic field extrapolation schemes. Due to the large area of missing data in
the example treated here in this work, the reconstructed field model obtained
must be treated with some caution. It is very likely that the field strength
in the area of missing data was small because the inversion procedure, which
calculates the surface field from the Stokes line spectra, abandons the
calculation if the signal is below a certain threshold. The magnetic field in
the corona, however, is dominated by the strongest flux elements on the
surface, even if they occupy only a small portion of the surface. We are
therefore confident that these dominant flux elements are accounted for in the
surface magnetogram, so that the resulting field model is fairly realistic. At
any rate, it is the field close to the real, which can be constructed from the
available sparse data. Therefore, we use a procedure which allows us to
incorporate measurement error and treat regions with lacking observational
data as in Tilaye:2010. The technique has been tested in cartesian geometry in
Wiegelmann10 for synthetic boundary data.

##  2 Optimization Principle in Spherical Geometry

Wheatland00 have proposed the variational principle to be solved iteratively
which minimizes Lorentz forces (one) and the divergence of magnetic field
(two) throughout the volume of interest, $V$. Later on the procedure has been
improved by Wiegelmann04 for cartesian geometry in a such way that it can only
uses the bottom boundary on the photosphere as an input. Here we use
optimization approach for the functional $(\mathcal{L}_{\mathrm{\omega}})$ in
spherical geometry [Wiegelmann (2007), Tadesse, Wiegelmann, and Inhester
(2009)] and iterate B to minimize $\mathcal{L}_{\mathrm{\omega}}$. The
modification concerns the input bottom boundary field
$\textbf{B}_{\textrm{obs}}$ which the model field B is not forced to match
exactly but we allow deviations of the order of the observational errors. The
modified variational problem is [Wiegelmann and Inhester (2010), Tadesse et
al. (2011)]:

$\textbf{B}=\textrm{argmin}(\mathcal{L}_{\omega})$
$\mathcal{L}_{\mathrm{\omega}}=\mathcal{L}_{\textrm{f}}+\mathcal{L}_{\textrm{d}}+\nu\mathcal{L}_{\textrm{photo}}\ilabel{four}$
(4)
$\mathcal{L}_{\textrm{f}}=\int_{V}\omega_{\textrm{f}}(r,\theta,\phi)B^{-2}\big{|}(\nabla\times\textbf{B})\times\textbf{B}\big{|}^{2}r^{2}\sin\theta
drd\theta d\phi$
$\mathcal{L}_{\textrm{d}}=\int_{V}\omega_{\textrm{d}}(r,\theta,\phi)\big{|}\nabla\cdot\textbf{B}\big{|}^{2}r^{2}\sin\theta
drd\theta d\phi$
$\mathcal{L}_{\textrm{photo}}=\int_{S}\big{(}\textbf{B}-\textbf{B}_{\textrm{obs}}\big{)}\cdot\textbf{W}(\theta,\phi)\cdot\big{(}\textbf{B}-\textbf{B}_{\textrm{obs}}\big{)}r^{2}\sin\theta
d\theta d\phi$

where $\mathcal{L}_{\mathrm{\textrm{f}}}$ and
$\mathcal{L}_{\mathrm{\textrm{d}}}$ measure how well the force-free Equations
(one) and divergence-free (two) conditions are fulfilled, respectively.
$\omega_{\textrm{f}}(r,\theta,\phi)$ and $\omega_{\textrm{d}}(r,\theta,\phi)$
are weighting functions for force-free term and divergence-free term,
respectively and identical for this study. The third integral,
$\mathcal{L}_{\textrm{photo}}$, is a surface integral over the photosphere
which serves to relax the field on the photosphere towards a force-free
solution without too much deviation from the original surface field data,
$\textbf{B}_{\textrm{obs}}$. In this integral,
$\textbf{W}(\theta,\phi)=\textrm{diag}(w_{\textrm{radial}},w_{\textrm{trans}},w_{\textrm{trans}})$
is diagonal matrix which gives different weights for observed surface field
components depending on its relative accuracy in measurement. In this sense,
lacking data is considered most inaccurate and is taken account of by setting
$W(\theta,\phi)$ to zero in all elements of the matrix.

Figure 1.: Surface contour plot of radial magnetic field component and vector
field plot of transverse field with white arrows. The color coding shows
$B_{r}$ on the photosphere. The vertical and horizontal axes show latitude,
$\theta$(in degree) and longitude, $\phi$(in degree) on the photosphere. In
the area coloured in olive, field values are lacking. The region inside the
black box corresponds to the physical domain where the weighting function is
unity and the outside region is the buffer zone where it declines to zero. The
blue boxes indicate the domains of the three active regions.fig1

We use a spherical grid $r$, $\theta$, $\phi$ with $n_{r}$, $n_{\theta}$,
$n_{\phi}$ grid points in the direction of radius, latitude, and longitude,
respectively. In the code, we normalize the magnetic field with the average
radial magnetic field on the photosphere and the length scale with a solar
radius for numerical reason. Figure fig1 shows a map of the radial component
of the field as color-coded and the transverse magnetic field depicted as
white arrows. For this particular dataset, about $86\%$ of the data pixels are
undetermined. The method works as follows:

  * •

We compute an initial source surface potential field in the computational
domain from $\textbf{B}_{\textrm{obs}}\cdot\hat{r}$, the normal component of
the surface field at the photosphere at $r=1R_{\mathrm{\odot}}$. The
computation is performed by assuming that a currentless ($\textbf{J}=0$ or
$\nabla\times\textbf{B}=0$) approximation holds between the photosphere and
some spherical surface $S_{s}$ (source surface where the magnetic field vector
is assumed radial). We compute the solution of this boundary-value problem in
a standard form of spherical harmonics expansion.

  * •

We minimize $\mathcal{L}_{\omega}$(Equations four) iteratively. The model
magnetic field B at the surface is gradually driven towards the observed field
$\textbf{B}_{\textrm{obs}}$ while the field in the volume $V$ relaxes to
force-free. If the observed field, $\textbf{B}_{\textrm{obs}}$, is
inconsistent, the difference $\textbf{B}-\textbf{B}_{\textrm{obs}}$ remains
finite depending in the control parameter $\nu$. At data gaps in
$\textbf{B}_{\textrm{obs}}$, the respective field value is automatically
ignored.

  * •

The iteration stops when $\mathcal{L}_{\omega}$ becomes stationary as
$\Delta\mathcal{L}_{\omega}/\mathcal{L}_{\omega}<10^{-4}$,
$\Delta\mathcal{L}_{\omega}$ is the decrease of $\mathcal{L}_{\omega}$ during
an iterative steps.

  * •

A convergence to $\mathcal{L}_{\omega}=0$ yields a perfect force-free and
divergence-free state and exact agreement of the boundary values B with
observations $\textbf{B}_{\textrm{obs}}$ in regions where the elements of W
are greater than zero. For inconsistent boundary data the force-free and
solenoidal conditions can still be fulfilled, but the surface term
$\mathcal{L}_{\textrm{photo}}$ will remain finite. This results in some
deviation of the bottom boundary data from the observations, especially in
regions where $w_{\textrm{radial}}$ and $w_{\textrm{trans}}$ are small. The
parameter $\nu$ is tuned so that these deviations do not exceed the local
estimated measurement error.

##  3 Results

In this work, we apply our extrapolation scheme to Milne-Eddington inverted
vector magnetograph data from the Synoptic Optical Long-term Investigations of
the Sun survey (SOLIS). As a first step for our work we remove non-magnetic
forces from the observed surface magnetic field using our spherical
preprocessing procedure. The code takes $\textbf{B}_{\textrm{obs}}$ as
improved boundary condition.

SOLIS/VSM provides full-disk vector-magnetograms, but for some individual
pixels the inversion from line profiles to field values may not have been
successful inverted and field data there will be missing for these pixels (see
Figure fig1). The different errors for the radial and transverse components of
$\textbf{B}_{\textrm{obs}}$ are taken account by different values for
$w_{\textrm{radial}}$ and $w_{\textrm{trans}}$. In this work we used
$w_{\textrm{radial}}=20w_{\textrm{trans}}$ for the surface preprocessed fields
as the radial component of $\textbf{B}_{\textrm{obs}}$ is measured with higher
accuracy.

We compute the 3D magnetic field above the observed surface region inside
wedge-shaped computational box of volume $V$, which includes an inner physical
domain $V^{\prime}$ and a buffer zone(the region outside the physical domain),
as shown in Figure fig1b. The physical domain $V^{\prime}$ is a wedge-shaped
volume, with two latitudinal boundaries at $\theta_{\textrm{min}}=-26^{\circ}$
and $\theta_{\textrm{max}}=16^{\circ}$, two longitudinal boundaries at
$\phi_{\textrm{min}}=129^{\circ}$ and $\phi_{\textrm{max}}=226^{\circ}$, and
two radial boundaries at the photosphere ($r=1R_{\odot}$) and
$r=1.75R_{\odot}$. Note that the longitude $\phi$ is measured from the center
meridian of the back side of the disk.

Figure 2.: Wedge-shaped computational box of volume $V$ with the inner
physical domain $V^{\prime}$ and a buffer zone. O is the center of the
Sun.fig1b

(a) 28 March 2008 15:45UT (b) 29 March 2008 15:48UT

(c) 30 March 2008 15:47UT (d) 28 March 2008 15:45UT

(e) 29 March 2008 15:48UT (f) 30 March 2008 15:47UT

(g) 28 March 2008 16:00UT (h) 29 March 2008 15:48UT

(i) 30 March 2008 15:48UT

Figure 3.: Top row: SOLIS/VSM magnetograms of respective dates. Middle row:
Magnetic field lines reconstructed from magnetograms on the top panel. Bottom
row: EIT image of the Sun at 195Å on indicated dates. fig2

(a) 28 March 2008 15:45UT (b) 29 March 2008 15:48UT (c) 30 March 2008 15:47UT

Figure 4.: Some magnetic field lines plots reconstructed from SOLIS
magnetograms using nonlinear force-free modelling. The color coding shows
$B_{r}$ on the photosphere.fig3

We define $V^{\prime}$ to be the inner region of $V$ (including the
photospheric boundary) with $\omega_{\textrm{f}}=\omega_{\textrm{d}}=1$
everywhere including its six inner boundaries $\delta V^{\prime}$. We use a
position-dependent weighting function to introduce a buffer boundary of
$nd=10$ grid points towards the side and top boundaries of the computational
box, $V$. The weighting functions, $\omega_{\textrm{f}}$ and
$\omega_{\textrm{d}}$ are chosen to be unity within the inner physical domain
$V^{\prime}$ and decline to 0 with a cosine profile in the buffer boundary
region [Wiegelmann (2004), Tadesse, Wiegelmann, and Inhester (2009)]. The
framed region in Figure fig1 corresponds to the lower boundary of the physical
domain $V^{\prime}$ with a resolution of $114\times 251$ pixels in the
photosphere.

The middle panel of Figure fig2 shows magnetic field line plots for three
consecutive dates of observation. The top and bottom panels of Figure fig2
show the position of the three active regions on the solar disk both for SOLIS
full-disk magnetogram 111http://solis.nso.edu/solis_data.html and SOHO/EIT
222http://sohowww.nascom.nasa.gov/data/archive image of the Sun observed at
195Å on the indicated dates and times. Figure fig3 shows some selected
magnetic field lines from reconstruction from the SOLIS magnetograms, zoomed
in from the middle panels of Figure fig2. In each column of Figure fig3 the
field lines are plotted from the same foot points to compare the change in
topology of the magnetic field over the period of the three days of
observation. In order to compare the fields at the three consecutive days
quantitatively, we computed the vector correlations between the three field
configurations. The vector correlation ($C_{\mathrm{\textrm{vec}}}$)
[Schrijver et al. (2006)] metric generalizes the standard correlation
coefficient for scalar functions and is given by

$C_{\mathrm{\textrm{vec}}}=\frac{\sum_{i}\textbf{v}_{i}\cdot\textbf{u}_{i}}{\sqrt{\sum_{i}|\textbf{v}_{i}|^{2}}\sqrt{\sum_{i}|\textbf{u}_{i}|^{2}}}\ilabel{nine}$
(5)

where $\textbf{v}_{i}$ and $\textbf{u}_{i}$ are 3D vectors at grid point $i$.
If the vector fields are identical, then $C_{\textrm{vec}}=1$; if
$\textbf{v}_{i}\perp\textbf{u}_{i}$ , then $C_{\textrm{vec}}=0$. The
correlations ($C_{\mathrm{\textrm{vec}}}$) of the 3D magnetic field vectors of
28 and 30 March with respect to the field on 29 March are $0.96$ and $0.93$
respectively. From these values we can see that there has been no major change
in the magnetic field configuration during this period.

(a) 28 March 2008 15:45UT (b) 29 March 2008 15:48UT (c) 30 March 2008 15:47UT

Figure 5.: Iso-surfaces (ISs) of the absolute NLFF magnetic energy density
($7.5\times 10^{16}$erg) for the three consecutive dates computed within the
entire computational domain.fig4

We also compute the values of the free magnetic energy estimated from the
excess energy of the extrapolated field beyond the potential field satisfying
the same $\textbf{B}_{\textrm{obs}}\cdot\hat{r}$ boundary condition. Similar
estimates have been made by Regnier and Thalmann for single active regions
observed at other times. From the corresponding potential and force-free
magnetic field, $\textbf{B}_{\textrm{pot}}$ and B, respectively, we can
estimate an upper limit to the free magnetic energy associated with coronal
currents

$E_{\mathrm{free}}=E_{\mathrm{\textrm{nlff}}}-E_{\mathrm{\textrm{pot}}}=\frac{1}{8\pi}\int_{V^{\prime}}\Big{(}B_{\textrm{nlff}}^{2}-B_{\textrm{pot}}^{2}\Big{)}r^{2}\textrm{sin}\theta
drd\theta d\phi.\ilabel{ten}$ (6)

Date | $E_{\textrm{nlff}}(10^{32}\textrm{erg})$ | $E_{\mathrm{\textrm{pot}}}(10^{32}\textrm{erg})$ | $E_{\mathrm{free}}(10^{32}\textrm{erg})$
---|---|---|---
28 March 2008 | $57.34$ | $53.89$ | $3.45$
29 March 2008 | $57.48$ | $54.07$ | $3.41$
30 March 2008 | $57.37$ | $53.93$ | $3.44$
Table 1.: The magnetic energy associated with extrapolated NLFF field
configurations for the three particular dates.table1

The computed energy values are listed in Table table1. The free energy on all
three days is about $3.5\times 10^{32}\textrm{erg}$. The magnetic energy
associated with the potential field configuration is about $54\times
10^{32}\textrm{erg}$. Hence $E_{\textrm{nlff}}$ exceeds $E_{\textrm{pot}}$ by
only 6$\%$. Figure fig4 shows Iso-surface plots of magnetic energy density in
the volume above the active regions. There are strong energy concentrations
above each active region. There were no major changes in the magnetic energy
density over the observation period and there was no major eruptive phenomenon
during those three days in the region observed.

In our previous work[Tadesse et al. (2011)], we have studied the connectivity
between two neighbouring active regions. In this work with an even larger
field of view, the three ARs share a decent amount of magnetic flux compared
to their internal flux from one polarity to the other (see Figure fig3). In
terms of the electric current they are much more isolated. In order to
quantify these connectivities, we have calculated the magnetic flux and the
electric currents shared between active regions. For the magnetic flux, e.g.,
we use

$\Phi_{\alpha\beta}=\sum_{i}|\textbf{B}_{i}\cdot\hat{r}|R^{2}_{\odot}\textrm{sin}(\theta_{i})\Delta\theta_{i}\Delta\phi_{i}\ilabel{ten}$
(7)

where the summation is over all pixels of $\mbox{AR}_{\alpha}$ from which the
field line ends in $\mbox{AR}_{\beta}$ or
$i\in\mbox{AR}_{\alpha}\|\,\mbox{conjugate footpoint}(i)\in\mbox{AR}_{\beta}$.
The indices $\alpha$ and $\beta$ denote the active regions and the index
number $1$ corresponds to AR $10989$, $2$ to AR $10988$ and $3$ to AR $10987$
of Figure fig1. For the electric current we replace the magnetic field,
$\textbf{B}_{i}\cdot\hat{r}$, by the vertical current density
$\textbf{J}_{i}\cdot\hat{r}$ in Equation (ten). Whenever the end point of a
field line falls outside (blue rectangles in Figure fig1) of the three ARs, we
categorize it as ending elsewhere. Both Table table2 and table3 show the
percentage of the total magnetic flux and electric current shared between the
three ARs. So, for example first column of Table table2 shows that $56.37\%$
of positive polarity of $\mbox{AR}_{1}$ is connected to negative polarity of
$\mbox{AR}_{1}$; line 2 shows that $13.66\%$ of positive/negative polarity of
$\mbox{AR}_{1}$ is connected to positive/negative polarity of $\mbox{AR}_{2}$,
and line 3 shows that there are no field lines $(0\%)$ connecting
positive/negative polarity of $\mbox{AR}_{1}$ with positive/negative polarity
of $\mbox{AR}_{3}$. The same technique applies for Table table3 too. The three
active regions are magnetically connected but much less by electric currents.

|  | $28th$ |  |  | $29th$ |  |  | $30th$ |
---|---|---|---|---|---|---|---|---|---
$\Phi_{\alpha\beta}$ | $\alpha=1$ | $2$ | $3$ | $\alpha=1$ | $2$ | $3$ | $\alpha=1$ | $2$ | $3$
$\beta=1$ | $56.37$ | $5.59$ | $0.00$ | $56.50$ | $5.48$ | $0.00$ | $56.50$ | $5.48$ | $0.00$
$2$ | $13.66$ | $81.12$ | $1.43$ | $13.66$ | $81.22$ | $1.43$ | $13.66$ | $81.22$ | $2.22$
$3$ | $0.00$ | $0.48$ | $71.47$ | $0.00$ | $0.48$ | $71.80$ | $0.00$ | $0.48$ | $71.80$
Elsewhere | $29.97$ | $12.82$ | $27.10$ | $29.84$ | $12.82$ | $26.77$ | $29.84$ | $12.82$ | $25.98$
Table 2.: The percentage of the total magnetic flux shared between the three ARs. $\Phi_{11}$, $\Phi_{22}$ and $\Phi_{33}$ denote magnetic flux of AR 10989(left), AR 10988(middle) and AR 10987(right) of Figure fig1, respectively. table2 |  | $28th$ |  |  | $29th$ |  |  | $30th$ |
---|---|---|---|---|---|---|---|---|---
$I_{\alpha\beta}$ | $\alpha=1$ | $2$ | $3$ | $\alpha=1$ | $2$ | $3$ | $\alpha=1$ | $2$ | $3$
$\beta=1$ | $82.47$ | $0.19$ | $0.00$ | $86.36$ | $0.19$ | $0.00$ | $94.16$ | $0.19$ | $0.00$
$2$ | $0.65$ | $85.25$ | $1.42$ | $0.65$ | $85.25$ | $1.42$ | $0.65$ | $85.25$ | $3.55$
$3$ | $0.00$ | $0.38$ | $82.27$ | $0.00$ | $0.38$ | $82.27$ | $0.00$ | $0.38$ | $82.27$
Elsewhere | $16.88$ | $14.18$ | $16.31$ | $12.99$ | $14.18$ | $16.31$ | $5.19$ | $14.18$ | $14.18$
Table 3.: The percentage of the total electric current shared between the
three ARs. $I_{11}$, $I_{22}$, and $I_{33}$ denote electric current of AR
10989(left), AR 10988(middle) and AR 10987(right) of Figure fig1,
respectively. table3

##  4 Conclusions

sect:disc We have investigated the coronal magnetic field associated with
three ARs 10987, 10987, 10989, on 28, 29 and 30 March 2008 by analysing
SOLIS/VSM data. We have used an optimization method for the reconstruction of
nonlinear force-free coronal magnetic fields in spherical geometry by
restricting the code to limited parts of the Sun [Wiegelmann (2007), Tadesse,
Wiegelmann, and Inhester (2009), Tadesse et al. (2011)]. The code was modified
so that it allows us to deal with lacking data and regions with poor signal-
to-noise ratio in a systematic manner[Wiegelmann and Inhester (2010), Tadesse
et al. (2011)].

We have studied the time evolution of magnetic field over the period of three
days and found no major changes in topologies as there was no major eruption
event. The magnetic energies calculated in the large wedge-shaped
computational box above the three active regions were not far apart in value.
This is the first study which contains three well separated ARs in our model.
This was made possible by the use of spherical coordinates and it allows us to
analyse linkage between the ARs. The active regions share a decent amount of
magnetic flux compared to their internal flux from one polarity to the other.
In terms of the electric current they are much more isolated.

## Acknowledgements

SOLIS/VSM vector magnetograms are produced cooperatively by NSF/NSO and
NASA/LWS. The National Solar Observatory (NSO) is operated by the Association
of Universities for Research in Astronomy, Inc., under cooperative agreement
with the National Science Foundation. Tilaye Tadesse Asfaw acknowledges a
fellowship of the International Max-Planck Research School at the Max-Planck
Institute for Solar System Research and the work of T. Wiegelmann was
supported by DLR-grant $50$ OC $0501$.

## References

  * Amari and Aly (2010) Amari, T., Aly, J.: 2010,. Astron. Astrophys. 522, A52.
  * Amari, Boulmezaoud, and Aly (2006) Amari, T., Boulmezaoud, T.Z., Aly, J.J.: 2006,. Astron. Astrophys. 446, 691\.
  * Amari, Boulmezaoud, and Mikic (1999) Amari, T., Boulmezaoud, T.Z., Mikic, Z.: 1999,. Astron. Astrophys. 350, 1051\.
  * Amari et al. (1997) Amari, T., Aly, J.J., Luciani, J.F., Boulmezaoud, T.Z., Mikic, Z.: 1997,. Solar Phys. 174, 129\.
  * Cuperman, Demoulin, and Semel (1991) Cuperman, S., Demoulin, P., Semel, M.: 1991,. Astron. Astrophys. 245, 285\.
  * Demoulin, Cuperman, and Semel (1992) Demoulin, P., Cuperman, S., Semel, M.: 1992,. Astron. Astrophys. 263, 351\.
  * DeRosa et al. (2009) DeRosa, M.L., Schrijver, C.J., Barnes, G., Leka, K.D., Lites, B.W., Aschwanden, M.J., et al.: 2009,. Astrophys. J. 696, 1780\.
  * Gary (2001) Gary, G.A.: 2001,. Solar Phys. 203, 71\.
  * Inhester and Wiegelmann (2006) Inhester, B., Wiegelmann, T.: 2006,. Solar Phys. 235, 201\.
  * Metcalf et al. (1995) Metcalf, T.R., Jiao, L., McClymont, A.N., Canfield, R.C., Uitenbroek, H.: 1995,. Astrophys. J. 439, 474\.
  * Metcalf et al. (2008) Metcalf, T.R., Derosa, M.L., Schrijver, C.J., Barnes, G., van Ballegooijen, A.A., Wiegelmann, T., Wheatland, M.S., Valori, G., McTtiernan, J.M.: 2008,. Solar Phys. 247, 269\.
  * Mikic and McClymont (1994) Mikic, Z., McClymont, A.N.: 1994,, Astronomical Society of the Pacific Conference Series 68, 225\.
  * Molodensky (1974) Molodensky, M.M.: 1974,. Solar Phys. 39, 393\.
  * Régnier and Priest (2007) Régnier, S., Priest, E.R.: 2007,. Astrophys. J. Lett. 669, L53.
  * Roumeliotis (1996) Roumeliotis, G.: 1996,. Astrophys. J. 473, 1095\.
  * Sakurai (1981) Sakurai, T.: 1981,. Solar Phys. 69, 343\.
  * Schrijver et al. (2006) Schrijver, C.J., Derosa, M.L., Metcalf, T.R., Liu, Y., McTiernan, J., Régnier, S., Valori, G., Wheatland, M.S., Wiegelmann, T.: 2006,. Solar Phys. 235, 161\.
  * Tadesse, Wiegelmann, and Inhester (2009) Tadesse, T., Wiegelmann, T., Inhester, B.: 2009,. Astron. Astrophys. 508, 421\.
  * Tadesse et al. (2011) Tadesse, T., Wiegelmann, T., Inhester, B., Pevtsov, A.: 2011,. Astron. Astrophys. 527, A30.
  * Thalmann, Wiegelmann, and Raouafi (2008) Thalmann, J.K., Wiegelmann, T., Raouafi, N.E.: 2008,. Astron. Astrophys. 488, L71.
  * Unno (1956) Unno, W.: 1956,. Publ. Astron. Soc. Japan 8, 108\.
  * Valori, Kliem, and Keppens (2005) Valori, G., Kliem, B., Keppens, R.: 2005,. Astron. Astrophys. 433, 335\.
  * Wheatland (2004) Wheatland, M.S.: 2004,. Solar Phys. 222, 247\.
  * Wheatland and Leka (2010) Wheatland, M.S., Leka, K.D.: 2010,. ArXiv e-prints.
  * Wheatland and Régnier (2009) Wheatland, M.S., Régnier, S.: 2009,. Astrophys. J. Lett. 700, L88.
  * Wheatland, Sturrock, and Roumeliotis (2000) Wheatland, M.S., Sturrock, P.A., Roumeliotis, G.: 2000,. Astrophys. J. 540, 1150\.
  * Wiegelmann (2004) Wiegelmann, T.: 2004,. Solar Phys. 219, 87\.
  * Wiegelmann (2007) Wiegelmann, T.: 2007,. Solar Phys. 240, 227\.
  * Wiegelmann (2008) Wiegelmann, T.: 2008,. J. Geophys. Res. 113, 3\.
  * Wiegelmann and Inhester (2010) Wiegelmann, T., Inhester, B.: 2010,. Astron. Astrophys. 516, A107.
  * Wiegelmann and Neukirch (2006) Wiegelmann, T., Neukirch, T.: 2006,. Astron. Astrophys. 457, 1053\.
  * Wiegelmann, Inhester, and Sakurai (2006) Wiegelmann, T., Inhester, B., Sakurai, T.: 2006,. Solar Phys. 233, 215\.
  * Wu et al. (1990) Wu, S.T., Sun, M.T., Chang, H.M., Hagyard, M.J., Gary, G.A.: 1990,. Astrophys. J. 362, 698\.
  * Yan and Sakurai (2000) Yan, Y., Sakurai, T.: 2000,. Solar Phys. 195, 89\.