# Leveraging Symbolic Algebra Systems to Simulate Contact Dynamics in Rigid
Body Systems

Simone Asci1 and Angadh Nanjangud2 *This work was not supported by any
organization1Simone Asci is with the School of Engineering and Materials
Science, Queen Mary University of London, Mile End Road, London, E1 4NS,
United Kingdom<EMAIL_ADDRESS>Nanjangud is with the School of
Engineering and Materials Science, Queen Mary University of London, Mile End
Road, London, E1 4NS, United Kingdom<EMAIL_ADDRESS>
###### Abstract

Collision detection plays a key role in the simulation of interacting rigid
bodies. However, owing to its computational complexity current methods
typically prioritize either maximizing processing speed or fidelity to real-
world behaviors. Fast real-time detection is achieved by simulating collisions
with simple geometric shapes whereas incorporating more realistic geometries
with multiple points of contact requires considerable computing power which
slows down collision detection. In this work, we present a new approach to
modeling and simulating collision-inclusive multibody dynamics by leveraging
computer algebra system (CAS). This approach offers flexibility in modeling a
diverse set of multibody systems applications ranging from human biomechanics
to space manipulators with docking interfaces, since the geometric
relationships between points and rigid bodies are handled in a generalizable
manner. We also analyze the performance of integrating this symbolic modeling
approach with collision detection formulated either as a traditional overlap
test or as a convex optimization problem. We compare these two collision
detection methods in different scenarios and collision resolution using a
penalty-based method to simulate dynamics. This work demonstrates an effective
simplification in solving collision dynamics problems using a symbolic
approach, especially for the algorithm based on convex optimization, which is
simpler to implement and, in complex collision scenarios, faster than the
overlap test.

##  I INTRODUCTION

Computer simulation of collision/contact dynamics is a germane research topic
of engineering science, particularly within the mechanics [1] and robotics [2]
communities. Current research on computational contact dynamics focuses on the
underlying numerical optimisation routines and addresses handling multiple
contacts for real-time simulation [3]. While important strides have been made
here, this research is limited in its employment of predefined models (e.g.,
6-DOF manipulators, humanoids, wheeled robots, quadrupeds), that do not
generalize to other domains such as spacecraft dynamics and control [4]. In
this paper, we present a new approach to simulating collision dynamics
multibody systems by integrating existing collision detection algorithms with
computer symbolic modeling for rigid body dynamics [5]; compared with
traditional approaches, the modeling capability facilitates the application of
the algorithm to uncommon and complicated shapes, enabling accurate results in
complex scenarios. Further, the symbolic approach is compatible with a variety
of contact dynamics models [6], which are appropriately formulated to describe
the system behavior through symbolic equations of motion. In our work, we
exploit an elastic-plastic contact model, commonly used in space manipulator
research [7], [8].

##  II SYMBOLIC APPROACH TO MODEL COLLISION

###  II-A Symbolic Simulation Framework

The proposed symbolic modeling approach utilizes SymPy, a widely used computer
algebra system (CAS) implemented in Python. We specifically make use of a
submodule that derives the symbolic equations of motion (EoMs) of multibody
systems [5]. The modular design of the framework consists of two parts: the
first models the dynamic system and generates the EoMs in symbolic form by
means of an automated routine; and the second piece converts symbolic EoMs to
their numerical equivalent that can be integrated over time to obtain the
evolution of the system.

###  II-B Collision Detection

Collision detection is defined as the procedure aimed to determine whether two
or more objects are overlapping. Specifically, in the context of dynamics
simulations, it detects when moving objects are in contact. It represents a
computational geometry problem with applications in various fields, including
computer graphics and video games. Popular algorithms for collision detection
are those based on Minkowsky difference, like the Gilbert-Johnson-Keerthi
(GJK) algorithm [9], [10]. Another category of collision detection algorithms
are based on the Separating Axis Theorem (SAT) [11], applied in synergy with
techniques to approximate the volume occupied by objects, like the Axis-
Aligned Bounding Box [12]. Such algorithms are employed in several physics
engines including Bullet [13], MuJoCo [14], and Box2D [15], where the trade-
off between detection accuracy and computational effort has a major impact on
the choice of the algorithm, since collision detection is mainly responsible
for simulation slowdown [4]. Specifically, the SAT is generally preferred for
simple applications, where accuracy plays a secondary role with respect to the
availability of computing power. Besides the aforementioned traditional
methods, examples of recently proposed methods are the Dynamic Collision
Checking which is based on heuristic search [16], and Fastron which leverages
machine learning techniques [17].

###  II-C Collision Resolution

Collision resolution is defined as the procedure to compute the dynamic
behavior of two or more bodies that are in a state of contact with each other.
This generally involves resolving the magnitude and direction of contact
forces and resulting moments exerted on the interacting bodies to then compute
the accelerations, though impulse-based approaches compute only the resulting
velocities of the bodies [18]. The collision force computation relates to
dynamics and depends primarily on the materials and mechanical properties of
the colliding bodies. Models proposed in the literature develop algebraic
models of collision for a pair of objects by running a series of collision
experiments and considering the relationship between the pre-impact and post-
impact velocities of each body[19]. These methods can be classified into two
categories: the discrete and continuous models [20]. The discrete approach,
known also as impulse-momentum or complementarity approach, is based on the
assumption that the poses of the bodies do not change significantly when
short-duration contact occurs; it models bodies as rigid and resolves contact
forces using kinematics constraints [21]. Complementarity methods handle non-
smooth events (e.g. collisions and contact interactions) by using impulsive
dynamics, unlike continuous methods that simulate body deformation during
contact and are affected by issues related to small time step size [22]. On
the other hand, the application of this approach to flexible and multibody
systems is complicated. The continuous approach, also referred to as force-
based approach or penalty-based method, approximates the local deformation of
the contacting bodies using the intersection between their respective
geometries, which is then utilised to model the contact force accordingly; it
has been widely applied in robotic contact problems, due to its suitability in
handling complex geometries [23], [24]. It includes the bristle-friction [25]
and the elastic-plastic deformation models [8]. The Bristle friction model is
based on a linear approximation of the Coulomb friction model. The elastic-
plastic model, which is utilized in this paper, calculates the interacting
force as the sum of the normal and tangent components to the impact surface,
the magnitudes of which depend on the relative velocity of the bodies and the
amount of local deformation.

###  II-D Collision response module architecture

Figure 1:  Architecture of the collision dynamics simulator

In a typical dynamics simulator, contact interactions are handled by a
collision response module (illustrated in Fig. 1 as the blue dashed box) which
takes as input the state of the bodies, i.e. their positions and velocities,
and returns a list of forces and moments to update the EoMs for each time
step. Its architecture is generally organized into a collision detection
module and a collision resolver module; the former detects any collision
between the bodies along with any other necessary information needed by the
resolver. This generally includes a metric that describes the distance between
the bodies, which is referred to as proximity (indicated in this paper with
parameter $\phi$, in SI units of $m$), and the minimum distance points
(referred to as MDPs in the paper), i.e. the points where the distance between
the two bodies is minimum and which therefore are potential future contact
points. In case of no overlap, the collision check is terminated and the
simulation continues without any update to the EoMs from the collision
response module. In this case $\phi$ still contains useful information that
can be utilized in the future time steps of the simulation. In case a
collision is detected, the collision detection module provides the required
information to the resolver which determines all the data necessary for the
simulator to resolve the outcome of the collision. Specifically, the
interpenetration (indicated with parameter $\rho$, in SI units of $m$) is the
metric equivalent to proximity and indicates the depth of penetration between
the bodies, while in this case the computed MDPs are actual contact points.
MDPs are used to compute the minimum translation vector (MTV) [26], which is
defined as the shortest distance along which objects should be moved away to
no longer be in a collision state. The direction of the contact force is
computed using the MTV, and its magnitude, which depends on the value of
$\rho$, is computed according to the specific contact method applied. Once the
interacting forces and moments are determined, the resolver module returns
them to the simulator which updates the states accordingly.

##  III Collision Detection Methods

###  III-A Separating Axis Theorem Based Collision Detection

Figure 2: Two convex shapes with their respective projections on a separating
axis

On a plane, the SAT states that two convex shapes do not intersect if and only
if there exists a line for which the projections of the objects onto this line
are not overlapping; this line is referred to as a separating axis (Fig. 2).
The same concept extends to the 3D case with a separating plane. The number of
lines to be tested depends on the shapes of the objects: two circles require
only a single test on the line joining the centers; while for polygons, the
normal to each edge of both shapes is a candidate separating axis. Projecting
a polygon onto a line requires performing the scalar product between each
vertex position and the unit vector lying along the line and storing the range
defined by minimum and maximum values; the same operation is repeated for the
vertices of the second object. The presence of an overlap can then be
determined by a comparison between the two ranges. Therefore the collision
detection check consists of a series of scalar products and can provide
information also on the amount of interpenetration between shapes, a useful
value to perform collision resolution at a later stage. All axes must be
tested for overlap to determine intersection, this aspect is responsible for
the algorithm’s computational complexity and significantly impacts
performance, especially when testing polygons with many edges or there are
many objects to test. However, at the first axis found where the projections
are not overlapping, the algorithm can immediately exit by concluding that the
shapes are not in a collision state. Therefore it is possible to speed up the
process by appropriately choosing the first axis to test [11], for example,
the line joining the centroids of the shapes has the highest probability of
detecting the absence of intersection between the shapes. The symbolic
approach is particularly advantageous for the SAT implementation because it
can handle symbolically the many geometric entities and linear algebra
operations required by this algorithm.

###  III-B Convex Optimization Based Collision Detection

Figure 3: Two convex shapes defined as the intersection of a set of half-
planes. The minimum distance points are highlighted in red.

Collision detection can also be formulated as an inequality-constrained convex
optimization (CO) program that performs a minimization over a continuous
objective function in the form [27]:

—l— x 12 x^TPx \+ c^Tx Gx ≤h

where $\mathbf{x}\in\mathbb{R}^{n}$ is the independent variable,
$\mathbf{P}\in\mathbb{S}_{+}^{n\times n}$ and $\mathbf{c}\in\mathbb{R}^{n}$
are quadratic and linear cost terms; the inequality constraint is described by
$\mathbf{G}\in\mathbb{R}^{l\times n}$, and $\mathbf{h}\in\mathbb{R}^{l}$. The
returned optimal value for the objective function indicates the proximity,
which is positive when no collision is detected and zero if the objects are
overlapping. The corresponding optimal values for the independent variables of
the problem indicate the coordinates of the MDPs (points $\tilde{A}$ and
$\tilde{B}$ in Fig. 3). The MDPs are constrained to be on the objects’
boundaries when there is no overlap, or in any point of the intersecting
region in case of collision. In this second case, indeed, any pair of
overlapping points located in aforementioned region satisfy the constraints of
belonging to both objects while having a distance equal to zero, which is an
optimal value for the objective function. The MDPs indicate the location where
the first contact occurred, consequently providing information to compute the
direction of collision force and the amount of interpenetration. Since, in the
case of objects intersection, the convex optimization algorithm returns
arbitrary values among a set of valid ones for the coordinates of the MDPs,
this means that convex optimization-based collision detection is not
compatible with penalty-based methods for collision resolution. Under certain
circumstances, it is possible to modify the problem to impose the position for
the MDPs on the boundaries and correctly compute collision resolution using a
method based on interpenetration. However, other methods could generally be
used, such as those based on non-penetration constraints, which model the
reaction force using Lagrangian multipliers [21]. In this work, cvxpy [28] is
utilized to solve the convex quadratic program as it offers a variety of
solvers based on methods commonly used, including the interior point and the
projected gradient methods. We make use of the embedded conic solver (ECOS),
an interior-point method for second-order cone programming.

##  IV Collision Between Rectangle and Circle

Figure 4:  Geometrical description of a rectangle and a circle as defined by
their centroids, orthogonal basis and parameters

This section describes a rectangle and circle in two-dimensional space (Fig.
4). To define the rectangle mathematically, let us consider a point $A^{*}$ in
the body-fixed reference frame
$\mathbf{A}=(\mathbf{\hat{a}_{1},\hat{a}_{2},\hat{a}_{3})}$ at a distance
$\mathbf{r_{A}}\in\mathbb{R}^{3}$ from $I^{*}$, the origin of an inertial
frame $\mathbf{I}=(\mathbf{\hat{i}_{1},\hat{i}_{2},\hat{i}_{3})}$. Note that
$\mathbf{\hat{i}_{3}}$ and $\mathbf{\hat{a}_{3}}$ are parallel and out-of-
plane. Then any point $y\in\mathbb{R}^{2}$ lying on the plane defined by
$\mathbf{\hat{a}_{1},\hat{a}_{2}}$ is within the rectangle if:

$\\{y\mid\mathbf{y}\cdot\mathbf{\hat{a}_{1}}\leq
C_{1}\wedge\mathbf{y}\cdot\mathbf{\hat{a}_{2}}\leq C_{2}\\}$ (1)

where $\mathbf{y}$ is the vector from origin $I^{*}$ to point $y$, parameters
$C_{1}$ and $C_{2}$ indicate rectangle dimensions, i.e. half-length and half-
width, and $\wedge$ is the logical conjunction operator.

Then, consider a circle with center $B^{*}$, of radius $R\in\mathbb{R}$
located at a distance $\mathbf{r_{B}}\in\mathbb{R}^{3}$ from origin $I^{*}$.
Its body-fixed reference frame is defined as
$\mathbf{B}=(\mathbf{\hat{b}_{1},\hat{b}_{2},\hat{b}_{3})}$, with
$\mathbf{\hat{b}_{3}}$ parallel to $\mathbf{\hat{i}_{3}}$. The circle is made
up of a set of points on the plane lying on axes $\mathbf{x_{I}}$ and
$\mathbf{y_{I}}$, such that any point $x\in\mathbb{R}^{3}$ is within a
distance $R$ from the center, as expressed in (2):

$\\{x\mid\left\|\mathbf{x}-\mathbf{r_{B}}\right\|\leq
R\wedge\mathbf{x}\cdot\mathbf{\hat{i}_{3}}=0\\}$ (2)

The following equation shows transformations between frames:

$\begin{split}&\mathbf{A}=\mathbf{{}^{I}{R}^{A}}\cdot\mathbf{I}\\\
&\mathbf{B}=\mathbf{{}^{I}{R}^{B}}\cdot\mathbf{I}\\\ \end{split}$ (3)

where $\mathbf{{}^{I}{R}^{A}}$ and $\mathbf{{}^{I}{R}^{B}}$ are the rotation
matrices relating the inertial frame $\mathbf{I}$ to frames $\mathbf{A}$ and
$\mathbf{B}$, respectively. They are defined as follows:

$\begin{split}&\mathbf{{}^{I}{R}^{A}}=\begin{bmatrix}\cos\theta_{1}&\sin\theta_{1}&0\\\
-\sin\theta_{1}&\cos\theta_{1}&0\\\ 0&0&1\end{bmatrix}\\\
&\mathbf{{}^{I}{R}^{B}}=\begin{bmatrix}\cos\theta_{2}&\sin\theta_{2}&0\\\
-\sin\theta_{2}&\cos\theta_{2}&0\\\ 0&0&1\end{bmatrix}\\\ \end{split}$ (4)

where $\theta_{1}$ and $\theta_{2}$ indicate the orientation angles. In the
remainder of this section, the vectors are expressed in the frame
$\mathbf{A}$, therefore the position of circle center $B^{*}$ is given by
$\mathbf{q}=\mathbf{{}^{I}{R}^{A}}(\mathbf{r_{B}-r_{A}})$.

###  IV-A Separating Axis Theorem Based Collision Detection

This procedure determines firstly the MDP of the rectangle, which is then used
to compute the MDP of the circle and the proximity value $\phi$. Let us
introduce two auxiliary quantities $\alpha$ and $\beta$:

$\begin{split}&\alpha=\lvert\mathbf{q}\cdot\mathbf{\hat{a}_{1}}\rvert-C_{1}\\\
&\beta=\lvert\mathbf{q}\cdot\mathbf{\hat{a}_{2}}\rvert-C_{2}\end{split}$ (5)

which indicate how far $B^{*}$ is from each rectangle side and in which region
of the space it is located; as shown in Fig. 5, point $B^{*}$ can be located
in:

  * •

corner regions and relative vertices, in green

  * •

upper and lower regions and relative edges, in blue

  * •

left and right regions and relative edges, in red

  * •

region inside the rectangle, in purple

Figure 5:  Rectangle regions

The MDP of the rectangle is indicated with $\tilde{A}$; it can be any point on
its boundary, including a vertex, depending on the circle center location. Its
coordinates $x$ and $y$ are defined by (6).

$\begin{cases}[\text{sgn}(\mathbf{q}{\cdot}\mathbf{\hat{a}_{1}})C_{1},\text{sgn}(\mathbf{q}{\cdot}\mathbf{\hat{a}_{2}})C_{2}]\
\text{for $(\alpha,\beta\geq 0)$}\\\
[\mathbf{q}{\cdot}\mathbf{\hat{a}_{1}},\text{sgn}(\mathbf{q}{\cdot}\mathbf{\hat{a}_{2}})C_{2}]\
\text{for $(\alpha<0\wedge\beta\geq 0)$}\\\
[\text{sgn}(\mathbf{q}{\cdot}\mathbf{\hat{a}_{1}})C_{1},\mathbf{q}{\cdot}\mathbf{\hat{a}_{2}}]\
\text{for $(\alpha\geq 0\wedge\beta<0)$}\\\
[\text{sgn}(\mathbf{q}{\cdot}\mathbf{\hat{a}_{1}})C_{1},\mathbf{q}{\cdot}\mathbf{\hat{a}_{2}}]\
\text{for $(\alpha,\beta<0\wedge\alpha>\beta)$}\\\
[\mathbf{q}{\cdot}\mathbf{\hat{a}_{1}},\text{sgn}(\mathbf{q}{\cdot}\mathbf{\hat{a}_{2}})C_{2}]\
\text{for $(\alpha,\beta<0\wedge\alpha<\beta)$}\\\
[\text{sgn}(\mathbf{q}{\cdot}\mathbf{\hat{a}_{1}})C_{1},\text{sgn}(\mathbf{q}{\cdot}\mathbf{\hat{a}_{2}})C_{2}]\
\text{for $(\alpha,\beta{<}0\wedge\alpha{=}\beta)$}\\\ \end{cases}$ (6)

Thus its position vector $\mathbf{\tilde{p}}$ is:

$\mathbf{\tilde{p}}=x\mathbf{\hat{a}_{1}}+y\mathbf{\hat{a}_{2}}$ (7)

Since the collision force is bound to point $\tilde{A}$, its position with
respect to the centroid of the rectangle is given by the vector
$\mathbf{p^{M}}$ (Fig. 4), which joins point $A^{*}$ to point $\tilde{A}$ and
coincides with $\mathbf{\tilde{p}}$. Consequently, the position of the MDP of
the circle $\tilde{B}$ with respect to $\mathbf{A}$ frame can be directly
computed as:

$\mathbf{\tilde{q}}=\mathbf{q}+R\frac{\mathbf{\tilde{p}}-\mathbf{q}}{\left\|\mathbf{\tilde{p}}-\mathbf{q}\right\|}$
(8)

The position $\mathbf{q^{M}}$ of the collision force applied to the circle
with respect to its center $B^{*}$ is:

$\mathbf{q^{M}}=R\frac{\mathbf{\tilde{p}}-\mathbf{q}}{\left\|\mathbf{\tilde{p}}-\mathbf{q}\right\|}$
(9)

and $\phi$ is defined as:

$\phi=\left\|\mathbf{\tilde{p}}-\mathbf{q}\right\|-R$ (10)

while $\rho$ is given by 11.

$\begin{cases}\rho=0&\text{for $\phi>0$}\\\ \rho=\lvert\phi\rvert&\text{for
$\phi\leq 0$}\end{cases}$ (11)

The components $n_{1},n_{2}$ of the normal to the surface can be calculated as
follows:

$\begin{cases}[\text{sgn}(\mathbf{q}{\cdot}\mathbf{\hat{a}_{1}})C_{1},\
\text{sgn}(\mathbf{q}{\cdot}\mathbf{\hat{a}_{2}})C_{2}]\ \text{for
$(\alpha,\beta\geq 0)$}\\\ [0,\
\text{sgn}(\mathbf{q}{\cdot}\mathbf{\hat{a}_{2}})]\ \text{for
$(\alpha<0\wedge\beta\geq 0)$}\\\
[\text{sgn}(\mathbf{q}{\cdot}\mathbf{\hat{a}_{1}}),\ 0]\ \text{for
$(\alpha\geq 0\wedge\beta<0)$}\\\
[\text{sgn}(\mathbf{q}{\cdot}\mathbf{\hat{a}_{1}}),\ 0]\ \text{for
$(\alpha,\beta<0\wedge\alpha>\beta)$}\\\ [0,\
\text{sgn}(\mathbf{q}{\cdot}\mathbf{\hat{a}_{2}})]\ \text{for
$(\alpha,\beta<0\wedge\alpha<\beta)$}\\\
[\text{sgn}(\mathbf{q}{\cdot}\mathbf{\hat{a}_{1}}),\
\text{sgn}(\mathbf{q}{\cdot}\mathbf{\hat{a}_{2}})]\ \text{for
$(\alpha,\beta<0\wedge\alpha=\beta)$}\\\ \end{cases}$ (12)

and finally the formulas of normal and tangent components of the force:

$\begin{split}&\mathbf{\hat{n}}=n_{1}\mathbf{\hat{a}_{1}}+n_{2}\mathbf{\hat{a}_{2}}\\\
&\mathbf{\hat{t}}=\mathbf{\hat{a}_{3}}\times\mathbf{\hat{n}}\end{split}$ (13)

###  IV-B Convex Optimization Based Collision Detection

Without loss of generality, all quantities are again expressed in the
rectangle reference frame. The optimization problem consists of finding two
points $\tilde{A}$ and $\tilde{B}$ belonging to the rectangle and circle,
respectively, minimizing the square of the distance between them. The
algorithm computes correctly the proximity and the MDPs coordinates only when
the objects do not intersect because in this case, every point inside the
overlapped regions is a candidate to be the MDP, the objective function is
always zero and therefore the value of the interpenetration cannot be
determined. However, a minor modification to the problem formulation makes it
possible to calculate $\rho$: by introducing a fictitious circle at the same
location as the original one, but with a smaller radius $R^{*}$ and solving
the following convex program for it:

—l— ~p, ~q^* ∥ ~p \- ~q^* ∥ ^2 ~p ⋅^a_1 ≤C_1 -~p ⋅^a_1 ≤C_1 ~p ⋅^a_2 ≤C_2 -~p
⋅^a_2 ≤C_2 ∥ q \- ~q^* ∥ ^2≤(R^*)^2

where $\mathbf{\tilde{q}^{*}}$ is the position of the MDP belonging to the
fictitious circle. Indicating with $b=R-R^{*}$ the difference between real and
fictitious radii, the interpenetration can be recovered from the surrogate
proximity $\phi^{*}$ with (14)

$\begin{cases}\rho=0&\text{for $\phi^{*}>b$}\\\ \rho=b-\phi^{*}&\text{for
$0<\phi^{*}\leq b$}\end{cases}$ (14)

$\phi^{*}$ is positive also when the original objects are in a collision
state, while it becomes zero if the rectangle and fictitious circle come into
contact, in which case the interpenetration can no longer be calculated
correctly. Vector $\mathbf{p^{M}}$ coincides with $\mathbf{\tilde{p}}$, which
is returned by the optimization algorithm; while vectors $\mathbf{\tilde{q}}$
and $\mathbf{q^{M}}$ are given by (8) and (9) respectively. Finally the normal
and tangent to the rectangle surface are computed as:

$\begin{split}&\mathbf{\hat{n}}=\frac{\mathbf{\tilde{q}^{*}}-\mathbf{\tilde{p}}}{\left\|\mathbf{\tilde{q}^{*}}-\mathbf{\tilde{p}}\right\|}\\\
&\mathbf{\hat{t}}=\mathbf{\hat{a}_{3}}\times\mathbf{\hat{n}}\end{split}$ (15)

##  V COLLISION RESOLUTION

In this study, collision resolution is implemented with a penalty-based
method. For the sake of simplicity, the reaction force applied at the
respective contact point is replaced by an equivalent force-moment system,
with the force applied at the center of mass of the body. The required
information provided by the collision detection algorithm is:

  * •

the amount of interpenetration

  * •

normal and tangent to the impact surface

  * •

the vector joining each body center of mass to the contact point

The magnitude of each force component is computed based on the elastic-plastic
approach [8], which defines the interaction between colliding bodies as a
spring-damped model. Complete information about material, geometry and
velocity of the bodies involved is supposed to be known. With reference to
(16) the normal force $F_{n}$ is composed of the elastic component, which is
proportional to the interpenetration value, and the plastic component, which
depends on the relative velocity along the normal direction. The tangent force
$F_{t}$ depends on the friction between the surfaces and the relative velocity
along the tangent direction.

$\begin{split}&F_{N}=(kc\cdot\rho^{3})\cdot(1-cc\cdot v_{N})\\\
&F_{T}=-mu\cdot F_{N}\cdot(\frac{2}{1+e^{\frac{-v_{T}}{vs}}}-1)\end{split}$
(16)

Where $kc$ and $cc$ are the contact stiffness parameter and damping
coefficient respectively, $mu$ is the friction coefficient and $vs$ is a
scaling factor, $v_{N}$ and $v_{T}$ are the normal and tangent values of the
relative velocity of the bodies. The computed force is applied, with opposite
signs, to both bodies; while the moment is obtained through a cross-product
with the respective force position vector.

##  VI RESULTS

Table I compares performances between the symbolic approach and the
traditional numeric one in simulating contact dynamics, applying SAT-based
method for collision detection and penalty-based method for collision
resolution. An additional comparison shows performances of the symbolic
approach integrated with CO-based collision detection. Five scenarios were
examined by running ten simulations each and calculating the average time, on
a computer equipped with an 8-core Intel Core i7-9750H CPU @ 2.60GHz.

TABLE I: SIMULATION TIME  Simulation scenario | Symbolic SAT | Numeric SAT | Symbolic CO
---|---|---|---
Bouncing Circle (2D) | 0.45 $\mathrm{s}$ | 0.53 $\mathrm{s}$ | 0.82 $\mathrm{s}$
Circle-Circle (2D) | 1.04 $\mathrm{s}$ | 0.93 $\mathrm{s}$ | 0.36 $\mathrm{s}$
Rectangle-Circle (2D) | 1.85 $\mathrm{s}$ | 1.40 $\mathrm{s}$ | 1.72 $\mathrm{s}$
Rectangle-Rectangle (2D) | 2.44 $\mathrm{s}$ | 2.54 $\mathrm{s}$ | 2.05 $\mathrm{s}$
Sphere-Cuboid (3D) | 4.01 $\mathrm{s}$ | 6.83 $\mathrm{s}$ | 3.09 $\mathrm{s}$

Regarding 2D scenarios, all methods achieve similar scores, but for the
numeric one, this comes at the cost of more modeling effort, which is not
reflected in the table. Convex optimization has a clear advantage in the 3D
case, where the SAT suffers from a significant increase in the number of tests
to be run. Furthermore, the implementation simplicity of CO algorithm scales
to more complex cases.

##  VII CONCLUSIONS

This paper presented a novel approach to implementing collision dynamics,
through the generation of symbolic equations of motion. The approach was
tested in five scenarios, integrating it with the Separating Axis Theorem and
convex optimization-based methods for collision detection of multibody systems
and with the elastic-plastic approach for collision resolution. We
demonstrated that it significantly simplifies the modeling process while
retaining the performance advantages of the standard numerical approach; this
allows modelers to investigate contact dynamics using a generalizable
approach. Further, a comparison between the two tested collision detection
methods highlighted that, for 2D cases, the SAT is fast and suitable for
applications with simple shapes and few contacts. On the other hand, the
convex optimization approach is better able to handle complex shapes, although
being generally not compatible with penalty-based methods. In applications
involving several contacts to check, this research indicates that convex
optimization-based collision detection is faster and more accurate than
algorithms based on the SAT.

## References

  * [1] T. A. Howell, S. L. Cleac’h, J. Z. Kolter, M. Schwager, and Z. Manchester, “Dojo: A differentiable simulator for robotics,” _arXiv preprint arXiv:2203.00806_ , 2022.
  * [2] H. Mylapilli and A. Jain, “Complementarity techniques for minimal coordinate contact dynamics,” _Journal of Computational and Nonlinear Dynamics_ , vol. 12, no. 2, 2017.
  * [3] M. Erdmann, “Multiple-point contact with friction: Computing forces and motions in configuration space,” in _Proceedings of 1993 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS’93)_ , vol. 1. IEEE, 1993, pp. 163–170.
  * [4] H. Choi, C. Crump, C. Duriez, A. Elmquist, G. Hager, D. Han, F. Hearl, J. Hodgins, A. Jain, F. Leve, _et al._ , “On the use of simulation in robotics: Opportunities, challenges, and suggestions for moving forward,” _Proceedings of the National Academy of Sciences_ , vol. 118, no. 1, 2021\.
  * [5] G. Gede, D. L. Peterson, A. S. Nanjangud, J. K. Moore, and M. Hubbard, “Constrained multibody dynamics with python: From symbolic equation generation to publication,” in _International Design Engineering Technical Conferences and Computers and Information in Engineering Conference_ , vol. 55973. American Society of Mechanical Engineers, 2013, p. V07BT10A051.
  * [6] E. Papadopoulos, F. Aghili, O. Ma, and R. Lampariello, “Robotic manipulation and capture in space: A survey,” _Frontiers in Robotics and AI_ , p. 228, 2021.
  * [7] H. Ghaednia, X. Wang, S. Saha, Y. Xu, A. Sharma, and R. L. Jackson, “A review of elastic–plastic contact mechanics,” _Applied Mechanics Reviews_ , vol. 69, no. 6, 2017.
  * [8] J. Oleś, T. Rybus, K. Seweryn, M. Surowiec, M. Wojtyra, M. Pietras, and M. Scheper, “Testing and simulation of contact during on-orbit operations,” in _Proceedings of 14th Symposium on Advanced Space Technologies in Robotics and Automation (ASTRA’2017)_. European Space Agency, 2017.
  * [9] E. G. Gilbert, D. W. Johnson, and S. S. Keerthi, “A fast procedure for computing the distance between complex objects in three-dimensional space,” _IEEE Journal on Robotics and Automation_ , vol. 4, no. 2, pp. 193–203, 1988\.
  * [10] “Enhancing gjk: Computing minimum and penetration distances between convex polyhedra,” in _Proceedings of international conference on robotics and automation_ , vol. 4. IEEE, 1997, pp. 3112–3117.
  * [11] S. Gottschalk, “Separating axis theorem. technical report tr96-024,” _Department of Computer Science. University of North Carolina, Chapel Hill_ , pp. 20–46, 1996.
  * [12] G. v. d. Bergen, “Efficient collision detection of complex deformable models using aabb trees,” _Journal of graphics tools_ , vol. 2, no. 4, pp. 1–13, 1997.
  * [13] M. Sagardia, T. Stouraitis, and J. Lopes e Silva, “A new fast and robust collision detection and force computation algorithm applied to the physics engine bullet: Method, integration, and evaluation,” in _Conference and Exhibition of the European Association of Virtual and Augmented Reality_ , 2014\.
  * [14] E. Todorov, T. Erez, and Y. Tassa, “Mujoco: A physics engine for model-based control,” in _2012 IEEE/RSJ international conference on intelligent robots and systems_. IEEE, 2012, pp. 5026–5033.
  * [15] I. Parberry, _Introduction to Game Physics with Box2D_. CRC Press, 2017.
  * [16] F. Schwarzer, M. Saha, and J.-C. Latombe, “Adaptive dynamic collision checking for single and multiple articulated robots in complex environments,” _IEEE Transactions on Robotics_ , vol. 21, no. 3, pp. 338–353, 2005.
  * [17] N. Das, N. Gupta, and M. Yip, “Fastron: An online learning-based model and active learning strategy for proxy collision detection,” in _Conference on Robot Learning_. PMLR, 2017, pp. 496–504.
  * [18] B. Mirtich and J. Canny, “Impulse-based simulation of rigid bodies,” in _Proceedings of the 1995 symposium on Interactive 3D graphics_ , 1995, pp. 181–ff.
  * [19] M. Mote, M. Egerstedt, E. Feron, A. Bylard, and M. Pavone, “Collision-inclusive trajectory optimization for free-flying spacecraft,” _Journal of Guidance, Control, and Dynamics_ , vol. 43, no. 7, pp. 1247–1258, 2020.
  * [20] G. Gilardi and I. Sharf, “Literature survey of contact dynamics modelling,” _Mechanism and machine theory_ , vol. 37, no. 10, pp. 1213–1239, 2002.
  * [21] D. E. Stewart and J. C. Trinkle, “An implicit time-stepping scheme for rigid body dynamics with inelastic collisions and coulomb friction,” _International Journal for Numerical Methods in Engineering_ , vol. 39, no. 15, pp. 2673–2691, 1996.
  * [22] A. Jain, “Minimal coordinates formulation of contact dynamics,” _Multibody dynamics_ , 2013.
  * [23] L. Luo and M. Nahon, “Development and validation of geometry-based compliant contact models,” _Journal of Computational and Nonlinear Dynamics_ , vol. 6, no. 1, 2011.
  * [24] S. Wu, F. Mou, and O. Ma, “Contact dynamics and control of a space manipulator capturing a rotating object,” in _AIAA Guidance, Navigation, and Control Conference_ , 2017, p. 1048.
  * [25] A. Peiret Giménez, J. Kövecses, and J. M. Font Llagunes, “The bristle friction model in a linear complementarity formulation for multibody systems with contacts,” in _Proceedings of the 4th Joint International Conference on Multibody System Dynamics_ , 2016, pp. 1–2.
  * [26] R. Mojtahedzadeh and A. J. Lilienthal, “A principle of minimum translation search approach for object pose refinement,” in _2015 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2015, pp. 2897–2903.
  * [27] K. Tracy, T. A. Howell, and Z. Manchester, “Diffpills: Differentiable collision detection for capsules and padded polygons,” _arXiv preprint arXiv:2207.00202_ , 2022.
  * [28] S. Diamond and S. Boyd, “CVXPY: A Python-embedded modeling language for convex optimization,” _Journal of Machine Learning Research_ , vol. 17, no. 83, pp. 1–5, 2016.


* (66) Alberto Martin-Lopez, Sergio Segura and Antonio Ruiz-Cortés  “A Catalogue of Inter-parameter Dependencies in RESTful Web APIs” Series Title: Lecture Notes in Computer Science  In _Service-Oriented Computing_ 11895 Cham: Springer International Publishing, 2019, pp. 399–414  DOI: 10.1007/978-3-030-33702-5˙31
  * (67) Alberto Martin-Lopez, Sergio Segura and Antonio Ruiz-Cortés  “Online testing of RESTful APIs: promises and challenges”  In _Proceedings of the 30th ACM Joint European Software Engineering Conference and Symposium on the Foundations of Software Engineering_ Singapore Singapore: ACM, 2022, pp. 408–420  DOI: 10.1145/3540250.3549144
  * (68) Alberto Martin-Lopez, Sergio Segura and Antonio Ruiz-Cortés  “RESTest: automated black-box testing of RESTful web APIs”  In _Proceedings of the 30th ACM SIGSOFT International Symposium on Software Testing and Analysis_ Virtual Denmark: ACM, 2021, pp. 682–685  DOI: 10.1145/3460319.3469082
  * (69) Alberto Martin-Lopez, Sergio Segura and Antonio Ruiz-Cortés  “RESTest: Black-Box Constraint-Based Testing of RESTful Web APIs” Series Title: Lecture Notes in Computer Science  In _Service-Oriented Computing_ 12571 Cham: Springer International Publishing, 2020, pp. 459–475  DOI: 10.1007/978-3-030-65310-1˙33
  * (70) Alberto Martin-Lopez, Sergio Segura and Antonio Ruiz-Cortés  “Test coverage criteria for RESTful web APIs”  In _Proceedings of the 10th ACM SIGSOFT International Workshop on Automating TEST Case Design, Selection, and Evaluation_ Tallinn Estonia: ACM, 2019, pp. 15–21  DOI: 10.1145/3340433.3342822
  * (71) Alberto Martin-Lopez, Sergio Segura, Carlos Müller and Antonio Ruiz-Cortés  “Specification and Automated Analysis of Inter-Parameter Dependencies in Web APIs” Conference Name: IEEE Transactions on Services Computing  In _IEEE Transactions on Services Computing_ 15.4, 2022, pp. 2342–2355  DOI: 10.1109/TSC.2021.3050610
  * (72) Barton P Miller, Lars Fredriksen and Bryan So  “An empirical study of the reliability of UNIX utilities”  In _Communications of the ACM_ 33.12 ACM New York, NY, USA, 1990, pp. 32–44
  * (73) A. Mirabella et al.  “Deep Learning-Based Prediction of Test Input Validity for RESTful APIs”  In _2021 IEEE/ACM Third International Workshop on Deep Learning for Testing and Testing for Deep Learning (DeepTest)_ Madrid, Spain: IEEE, 2021, pp. 9–16  DOI: 10.1109/DeepTest52559.2021.00008
  * (74) “OWASP API Security Project — OWASP Foundation”  URL: https://owasp.org/www-project-api-security/
  * (75) “OWASP API Security Project — OWASP Foundation”  URL: https://owasp.org/www-project-api-security/
  * (76) Lianglu Pan, Shaanan Cohney, Toby Murray and Van-Thuan Pham  “EDEFuzz: A Web API Fuzzer for Excessive Data Exposures”  In _Proceedings of the IEEE/ACM 46th International Conference on Software Engineering_ , ICSE ’24  New York, NY, USA: Association for Computing Machinery, 2024, pp. 1–12  DOI: 10.1145/3597503.3608133
  * (77) Chao Peng, Yujun Gao and Ping Yang  “Automated Server Testing: an Industrial Experience Report”  In _2022 IEEE International Conference on Software Maintenance and Evolution (ICSME)_ Limassol, Cyprus: IEEE, 2022, pp. 519–522  DOI: 10.1109/ICSME55016.2022.00071
  * (78) M. Pezzè and M. Young  “Software Testing and Analysis: Process, Principles, and Techniques”  Wiley India Pvt. Limited, 2008  URL: https://books.google.nl/books?id=7x10CgAAQBAJ
  * (79) L. Richardson, M. Amundsen and S. Ruby  “RESTful Web APIs: Services for a Changing World”  O’Reilly Media, 2013  URL: https://books.google.nl/books?id=ZXDGAAAAQBAJ
  * (80) José Miguel Rojas, Mattia Vivanti, Andrea Arcuri and Gordon Fraser  “A detailed investigation of the effectiveness of whole test suite generation”  In _Empirical Software Engineering_ 22.2, 2017, pp. 852–893  DOI: 10.1007/s10664-015-9424-2
  * (81) Orpheas Rooij et al.  “webFuzz: Grey-Box Fuzzing for Web Applications”  In _Computer Security – ESORICS 2021_ , Lecture Notes in Computer Science  Cham: Springer International Publishing, 2021, pp. 152–172  DOI: 10.1007/978-3-030-88418-5˙8
  * (82) Nicolás Serrano, Josune Hernantes and Gorka Gallardo  “Mobile Web Apps” Conference Name: IEEE Software  In _IEEE Software_ 30.5, 2013, pp. 22–27  DOI: 10.1109/MS.2013.111
  * (83) “Technical whitepaper for afl-fuzz”  In _Technical whitepaper for afl-fuzz_ URL: https://lcamtuf.coredump.cx/afl/technical_details.txt
  * (84) Erik Trickel et al.  “Toss a Fault to Your Witcher: Applying Grey-box Coverage-Guided Mutational Fuzzing to Detect SQL and Command Injection Vulnerabilities”  In _2023 IEEE Symposium on Security and Privacy (SP)_ , 44  SAN FRANCISCO: IEEE Computer Society, 2023, pp. 2658–2675  DOI: 10.1109/SP46215.2023.00007
  * (85) Omer Tripp et al.  “TAJ: Effective Taint Analysis of Web Applications”
  * (86) Chung-Hsuan Tsai, Shi-Chun Tsai and Shih-Kun Huang  “REST API Fuzzing by Coverage Level Guided Blackbox Testing” ISSN: 2693-9177  In _2021 IEEE 21st International Conference on Software Quality, Reliability and Security (QRS)_ IEEE, 2021, pp. 291–300  DOI: 10.1109/QRS54544.2021.00040
  * (87) Emanuele Viglianisi, Michael Dallago and Mariano Ceccato  “RESTTESTGEN: Automated Black-Box Testing of RESTful APIs” ISSN: 2159-4848  In _2020 IEEE 13th International Conference on Software Testing, Validation and Verification (ICST)_ IEEE, 2020, pp. 142–152  DOI: 10.1109/ICST46399.2020.00024
  * (88) “Vulnerability Scanning Tools — OWASP Foundation”  URL: https://owasp.org/www-community/Vulnerability_Scanning_Tools
  * (89) Khaing Hsu Wai et al.  “Code Modification Problems for Multimedia Use in JavaScript-Based Web Client Programming” Series Title: Lecture Notes in Networks and Systems  In _Complex, Intelligent and Software Intensive Systems_ 497 Cham: Springer International Publishing, 2022, pp. 548–556  DOI: 10.1007/978-3-031-08812-4˙53
  * (90) Seungwon Woo, Seungsoo Lee, Jinwoo Kim and Seungwon Shin  “RE-CHECKER: Towards Secure RESTful Service in Software-Defined Networking”  In _2018 IEEE Conference on Network Function Virtualization and Software Defined Networks (NFV-SDN)_ Verona, Italy: IEEE, 2018, pp. 1–5  DOI: 10.1109/NFV-SDN.2018.8725649
  * (91) Huayao Wu, Lixin Xu, Xintao Niu and Changhai Nie  “Combinatorial testing of RESTful APIs”  In _Proceedings of the 44th International Conference on Software Engineering_ Pittsburgh Pennsylvania: ACM, 2022, pp. 426–437  DOI: 10.1145/3510003.3510151
  * (92) Koji Yamamoto  “Efficient penetration of API sequences to test stateful RESTful services”  In _2021 IEEE International Conference on Web Services (ICWS)_ Chicago, IL, USA: IEEE, 2021, pp. 734–740  DOI: 10.1109/ICWS53863.2021.00101
  * (93) Rahulkrishna Yandrapally, Saurabh Sinha, Rachel Tzoref-Brill and Ali Mesbah  “Carving UI Tests to Generate API Tests and API Specification”  In _2023 IEEE/ACM 45th International Conference on Software Engineering (ICSE)_ Melbourne, Australia: IEEE, 2023, pp. 1971–1982  DOI: 10.1109/ICSE48619.2023.00167
  * (94) Joobeom Yun, Fayozbek Rustamov, Juhwan Kim and Youngjoo Shin  “Fuzzing of Embedded Systems: A Survey”  In _ACM Computing Surveys_ 55.7, 2023, pp. 1–33  DOI: 10.1145/3538644
  * (95) Michał Zalewski  “American Fuzzy Lop”, 2014  URL: http://lcamtuf.coredump.cx/afl
  * (96) Hao Zhang, Weiyu Dong and Liehui Jiang  “Zokfuzz: Detection of Web Vulnerabilities via Fuzzing”  In _2022 2nd International Conference on Consumer Electronics and Computer Engineering (ICCECE)_ IEEE, 2022, pp. 281–287  DOI: 10.1109/ICCECE54139.2022.9712748
  * (97) Jian Zhang et al.  “A Novel Neural Source Code Representation Based on Abstract Syntax Tree”  In _2019 IEEE/ACM 41st International Conference on Software Engineering (ICSE)_ Montreal, QC, Canada: IEEE, 2019, pp. 783–794  DOI: 10.1109/ICSE.2019.00086
  * (98) Man Zhang and Andrea Arcuri  “Adaptive Hypermutation for Search-Based System Test Generation: A Study on REST APIs with EvoMaster”  In _ACM Transactions on Software Engineering and Methodology_ 31.1, 2022, pp. 1–52  DOI: 10.1145/3464940
  * (99) Man Zhang and Andrea Arcuri  “Open Problems in Fuzzing RESTful APIs: A Comparison of Tools”  In _ACM Transactions on Software Engineering and Methodology_ 32.6, 2023, pp. 1–45  DOI: 10.1145/3597205
  * (100) Man Zhang, Asma Belhadi and Andrea Arcuri  “JavaScript SBST Heuristics to Enable Effective Fuzzing of NodeJS Web APIs”  In _ACM Transactions on Software Engineering and Methodology_ 32.6, 2023, pp. 1–29  DOI: 10.1145/3593801
  * (101) Man Zhang, Bogdan Marculescu and Andrea Arcuri  “Resource-based test case generation for RESTful web services”  In _Proceedings of the Genetic and Evolutionary Computation Conference_ Prague Czech Republic: ACM, 2019, pp. 1426–1434  DOI: 10.1145/3321707.3321815
  * (102) Man Zhang et al.  “White-Box Fuzzing RPC-Based APIs with EvoMaster: An Industrial Case Study”  In _ACM Transactions on Software Engineering and Methodology_ 32.5, 2023, pp. 122:1–122:38  DOI: 10.1145/3585009
  * (103) Jiazhen Zhao et al.  “Cefuzz: An Directed Fuzzing Framework for PHP RCE Vulnerability”  In _Electronics_ 11.5, 2022, pp. 758  DOI: 10.3390/electronics11050758
  * (104) Xiaogang Zhu, Sheng Wen, Seyit Camtepe and Yang Xiang  “Fuzzing: A Survey for Roadmap”  In _ACM Computing Surveys_ 54.11s, 2022, pp. 1–36  DOI: 10.1145/3512345

##  Appendix A OpenAPI Specification Example

⬇

1openapi: "3.1.0"

2info:

3 title: API Example

4paths:

5 /book:

6 get:

7 parameters:

8 - name: limit

9 in: query

10 description: How many items to return at one time (max 100)

11 required: false

12 schema:

13 type: integer

14 maximum: 100

15 format: int32

16 responses:

17 ’200’:

18 description: A paged array of books

19 headers:

20 x-next:

21 description: A link to the next page of responses

22 schema:

23 type: string

24 content:

25 application/json:

26 schema:

27 $ref: "#/components/schemas/books"

28 default:

29 description: unexpected error

30 content:

31 application/json:

32 schema:

33 $ref: "#/components/schemas/Error"

##  Appendix B Grammar example generated by RESTler atlidakis_restler_2019
(12)

⬇

15Request(

16 static("GET /api/customer/data"),

17 static("?idcustomer="),

18 fuzzable("integer"),

19 static("&name="),

20 fuzzable("string"),

21 static("HTTP/1.1"),

22 static("Accept:application/json"),

23 ......

24)


11institutetext: Department of Computer Science 22institutetext: Department of
Psychological & Brain Sciences
University of California, Santa Barbara, CA 93106 22email:
<EMAIL_ADDRESS>
# Greedy Optimization of Electrode Arrangement for Epiretinal Prostheses

Ashley Bruce 11 Michael Beyeler 1122 0000-0001-5233-844X

###### Abstract

Visual neuroprostheses are the only FDA-approved technology for the treatment
of retinal degenerative blindness. Although recent work has demonstrated a
systematic relationship between electrode location and the shape of the
elicited visual percept, this knowledge has yet to be incorporated into
retinal prosthesis design, where electrodes are typically arranged on either a
rectangular or hexagonal grid. Here we optimize the intraocular placement of
epiretinal electrodes using dictionary learning. Importantly, the optimization
process is informed by a previously established and psychophysically validated
model of simulated prosthetic vision. We systematically evaluate three
different electrode placement strategies across a wide range of possible
phosphene shapes and recommend electrode arrangements that maximize visual
subfield coverage. In the near future, our work may guide the prototyping of
next-generation neuroprostheses.

###### Keywords:

retinal prosthesis implant design dictionary selection.

##  1 Introduction

Current visual neuroprostheses consist of a microelectrode array (MEA)
implanted into the eye or brain that is used to electrically stimulate
surviving cells in the visual system in an effort to elicit visual percepts
(“phosphenes”). Current epiretinal implant users perceive highly distorted
percepts, which vary in shape not just across subjects, but also across
electrodes [7, 10], and may be caused by incidental stimulation of passing
nerve fiber bundles (NFBs) in the retina [4, 13]. However, this knowledge has
yet to be incorporated into prosthesis design.

To address this challenge, we make the following contributions:

  1. 1.

We optimize the intraocular placement of epiretinal electrodes using
dictionary learning. Importantly, this optimization process is informed by a
previously established and psychophysically validated model of simulated
prosthetic vision [1, 4].

  2. 2.

We systematically evaluate three different electrode placement strategies
across a wide range of possible phosphene shapes and recommend electrode
arrangements that maximize visual subfield coverage.

##  2 Related Work

Figure 1: A) Axon map model. _Left_ : Electrical stimulation (red disc) of a
NFB (gray lines) leads to tissue activation (dark-gray shaded region)
elongated along the NFB trajectory away from the optic disc (white circle).
The light-gray shaded region indicates the visual subfield that is being
simulated. _Right_ : The resulting visual percept appears elongated as well;
its shape can be described by two parameters, $\lambda$ (spatial extent along
the NFB trajectory) and $\rho$ (spatial extent perpendicular to the NFB). B)
As $\lambda$ increases, percepts become more elongated and start to overlap.

Sensory neuroprostheses such as retinal and cochlear implants are emerging as
a promising technology to restore lost sensory function. These devices bypass
the natural sensory transduction mechanism and provide direct electrical
stimulation of (retinal or auditory) nerve fibers that the brain interprets as
a (visual or auditory) percept. However, perceptual distortions can occur when
multiple electrodes stimulate the same neural pathways [2, 15], and electrode
placement has been shown to impact the quality of visual and hearing outcomes
[4, 8].

However, this information has yet to be incorporated into neuroprosthetic
device design and surgical placement, with most patients having less-than-
optimal MEA placement [3, 6]. Whereas previous studies have optimized the
shape of individual electrodes [12], which electrodes to activate in order to
produce a desired visual response [5, 9, 14], or the overall implant placement
[3], we are unaware of any studies that have attempted to optimize the
placement of individual electrodes within an implant.

In the case of retinal implants, recent work has demonstrated that phosphene
shape strongly depends on the retinal location of the stimulating electrode
[4]. Because retinal ganglion cells (RGCs) send their axons on highly
stereotyped pathways to the optic nerve, an electrode that stimulates nearby
axonal fibers would be expected to antidromically activate RGC bodies located
peripheral to the point of stimulation, leading to percepts that appear
elongated in the direction of the underlying NFB trajectory (Fig. 1A,
_right_). Using a simulated map of NFBs, Reference [4] was thus able to
accurately predict phosphene shape for various users of the Argus Retinal
Prosthesis System (Second Sight Medical Products, Inc.), by assuming that an
axon’s sensitivity to electrical stimulation:

  1. i.

decays exponentially with decay constant $\rho$ as a function of distance from
the stimulation site,

  2. ii.

decays exponentially with decay constant $\lambda$ as a function of distance
from the cell body, measured as axon path length.

As can be seen in Fig. 1B, electrodes near the horizontal meridian are
predicted to elicit circular percepts, while other electrodes are predicted to
produce elongated percepts that will differ in angle based on whether they
fall above or below the horizontal meridian. In addition, the values of $\rho$
and $\lambda$ dictate the size and elongation of elicited phosphenes,
respectively, which may drastically affect visual outcomes. Specifically, if
two electrodes happen to activate the same NFB, they might not generate two
distinct phosphenes (Fig. 1B, _right_).

Instead of arranging electrodes such that they efficiently tile the _retinal
surface_ , a better approach might thus be to arrange electrodes such that the
_elicited percepts_ effectively tile the _visual field_. In the following, we
will demonstrate that this is equivalent to finding the smallest set of
electrodes that cover a desired portion of the visual field (here termed a
visual subfield). However, this strategy can be readily applied wherever there
is a topological mapping from stimulus space to perceptual space (_e.g._ ,
visual, auditory, tactile stimulation). Incorporating this knowledge into
implant design could therefore be indispensable to the success of future
visual neuroprostheses.

##  3 Methods

###  3.1 Phosphene Model

Let $\mathcal{E}=\\{e_{1},\ldots,e_{N}\\}$ be the set of $N$ electrodes in a
MEA, where the $i$-th electrode $e_{i}=(x_{i},y_{i},r_{i})$ is described by
its location on the retinal surface $(x_{i},y_{i})$ and its radius $r_{i}>0$.
For example, Argus II can be described by $\mathcal{E}_{\mathrm{ArgusII}}$,
where $N=|\mathcal{E_{\mathrm{ArgusII}}}|=60$,
$r_{i}=$122.5\text{\,}\mathrm{\SIUnitSymbolMicro m}$$ $\forall i$, and
$(x_{i},y_{i})$ are spaced $575\text{\,}\mathrm{\SIUnitSymbolMicro m}$ apart
on a rectangular grid. We do not assume any particular ordering of
$\mathcal{E}$.

Furthermore, let $\mathcal{S}=\\{(s_{1},\ldots,s_{k})\\}$ be the set of
stimuli sent to $k\leq N$ electrodes in the MEA, where the $i$-th stimulus
$s_{i}=(e_{i},a_{i})$ is described by an electrode $e_{i}\in\mathcal{E}$ and
its corresponding activation function $a_{i}$. In practice, $a_{i}$ may be a
biphasic pulse train of a given duration, pulse amplitude, and pulse
frequency. However, for the purpose of this study, we limited ourselves to
spatial activation values $a_{i}\in\mathbb{R}_{\geq 0}$, which did not contain
a temporal component.

A phosphene model $\mathcal{M}$ then takes a set of stimuli $\mathcal{S}$ as
input and outputs a visual percept $p\in\mathbb{R}^{H\times W}_{\geq 0}$,
which is a height ($H$) $\times$ width ($W$) grayscale image. In general,
$\mathcal{M}$ is a nonlinear function of $\mathcal{S}$ and depends on subject-
specific parameters $\theta_{\mathrm{subject}}$ such as $\rho$ and $\lambda$;
thus $p=\mathcal{M}(\mathcal{S};\theta_{\mathrm{subject}})$.

For the purpose of this study, we used _pulse2percept_ 0.8.0, a Python-based
simulation framework for bionic vision [1] that provides an open-source
implementation of the axon map model [4], as described in the previous
section. Constrained by electrophysiological and psychophysical data, this
model predicts what a bionic eye user should “see” for any given set of
stimuli $\mathcal{S}$.

###  3.2 Dictionary Selection

####  3.2.1 Problem Formulation

Let $\mathcal{T}$ now be the set of all nonoverlapping epiretinal electrodes,
and $\mathcal{D}$ be a subset of those; that is,
$\mathcal{D}\subset\mathcal{T}$, where $|\mathcal{D}|\ll|\mathcal{T}|$. In the
dictionary selection problem, we are interested in finding the dictionary
$\mathcal{D}^{*}$ that maximizes a utility function $F$ (see next subsection):

$\mathcal{D^{*}}=\operatorname*{arg\,max}_{|\mathcal{D}|\leq
k}F(\mathcal{D}),$ (1)

where $k$ is a constraint on the number of electrodes that the dictionary can
be composed of.

This optimization problem presents combinatorial challenges, as we have to
find the set $\mathcal{D}^{*}$ out of exponentially many options in
$\mathcal{T}$. However, we will only consider utility functions $F$ with the
following properties:

  1. i.

The empty set has zero utility; that is, $F(\emptyset)=0$.

  2. ii.

$F$ increases monotonically; that is, whenever
$\mathcal{D}\subseteq\mathcal{D^{\prime}}$, then $F(\mathcal{D})\leq
F(\mathcal{D^{\prime}})$.

  3. iii.

$F$ is approximately submodular; that is, there exists an $\varepsilon$ such
that whenever $\mathcal{D}\subseteq\mathcal{D^{\prime}}\subseteq\mathcal{T}$
and an electrode $e\in\mathcal{T}\setminus\mathcal{D^{\prime}}$, it holds that
$F(\mathcal{D}\cup\\{e\\})-F(\mathcal{D})\geq
F(\mathcal{D^{\prime}}\cup\\{e\\})-F(\mathcal{D^{\prime}})-\varepsilon$. This
property implies that adding a new electrode $e$ to a larger dictionary
$\mathcal{D^{\prime}}$ helps at most $\varepsilon$ more than adding $e$ to a
subset $\mathcal{D}\subseteq\mathcal{D^{\prime}}$.

For utility functions with the above properties, Nemhauser et al. [11] proved
that a simple greedy algorithm that starts with the empty set
$\mathcal{D}_{0}=\emptyset$, and at every iteration $i$ adds the element

$d_{i}=\operatorname*{arg\,max}_{d\in\mathcal{T}\setminus\mathcal{D}}F(\mathcal{D}_{i-1}\cup\\{d\\}),$
(2)

where $\mathcal{D}_{i}=\\{d_{1},\ldots,d_{i}\\}$, is able to obtain a near-
optimal solution.

####  3.2.2 Utility Function

Ideally, the utility function $F$ would directly assess the quality of the
generated artificial vision. As a first step towards such a quality measure,
we considered the ability of a set of electrodes to lead to phosphenes that
cover a specific visual subfield (i.e., the gray shaded region in Fig. 1). We
would thus activate every electrode in $\mathcal{D}$, represented by
$\mathcal{S}_{\mathcal{D}}$, and calculate the percept
$p_{\mathcal{D}}=\mathcal{M}(\mathcal{S}_{\mathcal{D}};\theta_{\mathrm{subject}})$.
Then $F$ was given as the visual subfield coverage; that is:

$F(\mathcal{D})=\sum_{w=1}^{w=W}\sum_{h=1}^{h=H}p_{hw}\geq\varepsilon,$ (3)

where $W$ and $H$ were the width and height of the percept, respectively,
$\varepsilon=0.1\max_{\forall h,w}(p_{hw}$), and $F\in[0,HW]$.

####  3.2.3 Dictionary Selection Strategies

To find $\mathcal{D}^{*}$, we considered three different strategies:

  * •

No Overlap: At each iteration $i$, $d_{i}$ was chosen according to Eq. 2.

  * •

Not Too Close: To consider manufacturing constraints, we modified Eq. 2 above
to enforce that all electrodes were placed at least
$c=$112.5\text{\,}\mathrm{\SIUnitSymbolMicro m}$$ apart from each other:

$d_{i}=\operatorname*{arg\,max}_{d\in\mathcal{T}\setminus\mathcal{D}\;\mathrm{s.t.}\;||d,d_{j}||_{2}\geq
c\;\forall d_{j}\>\in\>\mathcal{D}\setminus
d}F(\mathcal{D}_{i-1}\cup\\{d\\}).$ (4)

  * •

Max Pairwise Distance: As electrical crosstalk is one of the main causes of
impaired spatial resolution in retinal implants [15], we further modified Eq.
4 to place electrodes as far away from each other as possible:

$d_{i}=\operatorname*{arg\,max}_{d\in\mathcal{T}\setminus\mathcal{D}\;\mathrm{s.t.}\;||d,d_{j}||_{2}\geq
c\;\forall d_{j}\>\in\>\mathcal{D}\setminus
d}F(\mathcal{D}_{i-1}\cup\\{d\\})+\alpha\sum_{d_{i}\in\mathcal{D}\setminus
d}||d,d_{i}||_{2},$ (5)

where $\alpha=$1\text{\times}{10}^{-4}$$ was a scaling factor.

####  3.2.4 Implementation Details

For the sake of feasibility, we limited $\mathcal{T}$ to electrodes placed on
a finely spaced search grid,
$x_{i}\in[-3000,3000]\>$\text{\,}\mathrm{\SIUnitSymbolMicro m}$$ and
$y_{i}\in[-2000,2000]\>$\text{\,}\mathrm{\SIUnitSymbolMicro m}$$, sampled at
$112.5\text{\,}\mathrm{\SIUnitSymbolMicro m}$ (i.e., the radius of an Argus II
electrode). This led to a manageable set size ($|\mathcal{T}|\approx 2000$)
while still allowing electrodes to be placed right next to each other (if
desirable). To find the electrode $d_{i}$ at each iteration in Eqs. 2, 4, and
5 above, we thus performed a grid search.

####  3.2.5 Stopping Criteria

The dictionary search was stopped when at least one of the following criteria
were met:

  * •

visual subfield coverage reached $99\text{\,}\mathrm{\char 37\relax}$,

  * •

the utility score $F$ did not improve by $\geq$1\text{\times}{10}^{-6}$$ on
two consecutive runs,

  * •

no more viable electrode locations were available (i.e.,
$\mathcal{T}\setminus\mathcal{D}\;\mathrm{s.t.}\;||d,d_{i}||_{2}\geq
c\;\forall d_{i}\>\in\>\mathcal{D}\setminus d=\emptyset$).

##  4 Results

###  4.1 Visual Subfield Coverage

The results of the greedy dictionary selection are shown in Fig. 2. For all
three dictionary selection strategies, the number of electrodes required to
cover at least $99\text{\,}\mathrm{\char 37\relax}$ of the visual subfield was
inversely proportional to $\rho$ and $\lambda$. As expected, the largest
number was achieved with the smallest, most compact phosphene shape
($\rho=$100\text{\,}\mathrm{\SIUnitSymbolMicro
m}$,\lambda=$200\text{\,}\mathrm{\SIUnitSymbolMicro m}$$). The required
electrode number dropped rapidly with increasing $\rho$ and $\lambda$,
indicating that for large phosphenes, a prototype implant such as Argus I
($4\times 4$ electrodes) might be sufficient to cover the whole subfield.

It is interesting to note that, for any given $\rho$ and $\lambda$
combination, the Max Pairwise Distance strategy required a smaller number of
electrodes than the No Overlap strategy. However, using the Not Too Close
strategy, smaller $\rho$ and $\lambda$ combinations were no longer able to
reach full coverage, as electrodes could no longer be placed too close to each
other or to the visual subfield boundary. This issue was amplified with the
Max Pairwise Distance strategy, with which coverage dropped for most $\rho$
and $\lambda$ combinations to $95\text{\,}\mathrm{\char 37\relax}$.

Figure 2: Number of electrodes needed to reach maximum visual subfield
coverage for the three dictionary selection strategies ($\dagger$:
$80-$95\text{\,}\mathrm{\char 37\relax}$$ coverage, $*$:
$95-$99\text{\,}\mathrm{\char 37\relax}$$ coverage). Max Pairwise Distance
always had $*$, unless otherwise noted by $\dagger$. Figure 3: Representative
examples of the final electrode arrangements generated with the three
dictionary selection strategies for different $\rho$ and $\lambda$
combinations. The left half of each panel shows the retinal location of all
electrodes (small circles) in the implant, and the shaded region indicates the
visual subfield (compare to Fig. 1). The visual percept that results from
simultaneously activating every electrode in the implant is shown in the right
half of each panel.

With larger $\rho$ and $\lambda$ values, an implant required less than ten
electrodes to cover the visual subfield. However, such a small number was not
necessarily desirable, as it also reduced the number of distinct phosphenes
that the implant can produce. Rather than focusing on visual subfield coverage
alone, one might therefore ask what kinds of electrode arrangements the three
dictionary selection strategies yield and what the resulting percepts look
like.

###  4.2 Electrode Arrangement

Example electrode arrangements suggested by the three dictionary selection
strategies are shown in Fig. 3. Here it is evident that, as $\lambda$
increased, electrodes were preferentially placed on the top, bottom, and right
boundaries of the visual subfield. This placement would often lead to the
longest streaks, thus yielding the largest coverage. For small $\rho$ values,
electrodes aggregated mainly on NFBs, spaced $\lambda$ apart, so that the
streaks generated by different electrodes tiled the visual subfield. As $\rho$
increased, electrodes tended to migrate away from the border and more inward,
due to the outward spread of the generated percept.

Figure 4: A) Representative samples from the MNIST database (_top row_)
represented with an optimized 60-electrode implant (“No Overlap”, _middle
row_) and Argus II (_bottom row_). B) Mean absolute error for 100 randomly
selected MNIST digits (vertical bars: standard error of the mean, ***:
$p<.001$). C) Number of active electrodes needed to represent the 100 digits.
All simulations had $\rho=$200\text{\,}\mathrm{\SIUnitSymbolMicro
m}$,\lambda=$400\text{\,}\mathrm{\SIUnitSymbolMicro m}$$.

###  4.3 Comparison with Argus II

To assess whether the optimized electrode arrangements provided an improvement
over a rectangular MEA, we compared our results to Argus II. We thus modified
our experiment such that the dictionary consisted of at most 60 electrodes
(i.e., same as Argus II). We then considered the ability of the optimized
implant to represent handwritten digits from the MNIST database (Fig. 4A). We
found that the optimized implant not only led to smaller errors between
predicted percepts and ground-truth digits (paired $t$-test, $p<.001$; Fig.
4B) but also required less active electrodes overall (paired $t$-test,
$p<.001$; Fig. 4C).

Furthermore, we calculated the visual subfield coverage achieved with these 60
electrodes and compared the result to Argus II (Table. 1). Given the fixed
number of electrodes, our electrode selection strategies were always able to
cover a larger portion of the visual subfield than Argus II.

Overall these results suggest that a rectangular grid might not be the best
electrode arrangement for epiretinal prostheses.

##  5 Conclusion

We report epiretinal electrode arrangements that maximize visual subfield
coverage as discovered by dictionary learning. We were able to obtain nearly
full coverage across different $\rho$ and $\lambda$ values for a given implant
size and electrode radius, and report the number of electrodes necessary to
reach this coverage. Future work should extend the dictionary selection
strategy to other implant sizes and explore all possible stimuli, not just the
ones that may result in sufficient visual subfield coverage.

This preliminary study is a first step towards the use of computer simulations
in the prototyping of novel neuroprostheses. To the best of our knowledge,
this is the first study using a psychophysically validated phosphene model to
optimize electrode arrangement. Even though we have focused on a specific
implant technology, our strategy can be readily applied wherever there is a
topological mapping from stimulus space to perceptual space (e.g., visual,
auditory, tactile). This means that our approach could be extended to other
neuromodulation technologies that include (but are not limited to) other
electronic prostheses and optogenetic technologies. In the near future, our
work may therefore guide the prototyping of next-generation epiretinal
prostheses.

$\rho$ | $100\text{\,}\mathrm{\SIUnitSymbolMicro m}$ | $150\text{\,}\mathrm{\SIUnitSymbolMicro m}$ | $200\text{\,}\mathrm{\SIUnitSymbolMicro m}$
---|---|---|---
$\lambda$ | AII | NTC | MPD | AII | NTC | MPD | AII | NTC | MPD
$200\text{\,}\mathrm{\SIUnitSymbolMicro m}$ | 39.2 | 43.2 | 42.9 | 74.1 | 75.6 | 75.9 | 86.2 |  $>$99* | 95.2
$400\text{\,}\mathrm{\SIUnitSymbolMicro m}$ | 55.0 | 66.8 | 66.8 | 80.6 | 91.3 | 91.4 | 87.5 |  $>$99* | 95.1
$600\text{\,}\mathrm{\SIUnitSymbolMicro m}$ | 62.4 | 81.8 | 81.3 | 81.9 | 96.9 | 95.3 | 88.4 |  $>$99* | 95.2
$800\text{\,}\mathrm{\SIUnitSymbolMicro m}$ | 66.7 | 88.8 | 88.8 | 82.8 | 97.7 | 95.2 | 89.0 |  $>$99* | 95.4
$1000\text{\,}\mathrm{\SIUnitSymbolMicro m}$ | 69.4 | 92.5 | 92.2 | 83.5 | 98.1* | 95.2 | 89.6 |  $>$99* | 95.8
$1200\text{\,}\mathrm{\SIUnitSymbolMicro m}$ | 71.9 | 92.6 | 92.5 | 84.2 | 98.0* | 95.3 | 90.1 |  $>$99* | 95.7
$1400\text{\,}\mathrm{\SIUnitSymbolMicro m}$ | 72.4 | 93.9* | 93.8 | 84.6 | 98.1* | 95.27 | 90.4 |  $>$99* | 95.6
$1600\text{\,}\mathrm{\SIUnitSymbolMicro m}$ | 73.4 | 93.9* | 93.8 | 85.0 | 98.5* | 95.41 | 90.7 |  $>$99* | 95.8
$1800\text{\,}\mathrm{\SIUnitSymbolMicro m}$ | 74.1 | 94.0* | 93.9 | 85.4 | 98.6* | 95.16 | 91.0 |  $>$99* | 95.7
$2000\text{\,}\mathrm{\SIUnitSymbolMicro m}$ | 74.7 | 94.4* | 94.4 | 85.9 | 98.3* | 95.08 | 91.4 |  $>$99* | 95.2

Table 1: Visual subfield coverage achieved with 60 electrodes using different
electrode placement strategies. AII: Argus II. NTC: Not Too Close. MPD: Max
Pairwise Distance. Asterisk (*): Search terminated before utilizing 60
electrodes.

## Acknowledgements

This work was supported by the National Institutes of Health (NIH R00
EY-029329 to MB). The authors would like to thank Madori Spiker and Alex Rasla
for valuable contributions to an earlier version of the project.

## References

  * [1] Beyeler, M., Boynton, G.M., Fine, I., Rokem, A.: pulse2percept: A Python-based simulation framework for bionic vision. In: Huff, K., Lippa, D., Niederhut, D., Pacer, M. (eds.) Proceedings of the 16th Science in Python Conference. pp. 81–88 (2017), https://doi.org/10.25080/shinma-7f4c6e7-00c
  * [2] Beyeler, M., Rokem, A., Boynton, G.M., Fine, I.: Learning to see again: biological constraints on cortical plasticity and the implications for sight restoration technologies. J Neural Eng 14(5), 051003 (Jun 2017), https://doi.org/10.1088/1741-2552/aa795e
  * [3] Beyeler, M., Boynton, G.M., Fine, I., Rokem, A.: Model-Based Recommendations for Optimal Surgical Placement of Epiretinal Implants. In: Shen, D., Liu, T., Peters, T.M., Staib, L.H., Essert, C., Zhou, S., Yap, P.T., Khan, A. (eds.) Medical Image Computing and Computer Assisted Intervention – MICCAI 2019. pp. 394–402. Lecture Notes in Computer Science, Springer International Publishing (2019), https://doi.org/10.1007/978-3-030-32254-0˙44
  * [4] Beyeler, M., Nanduri, D., Weiland, J.D., Rokem, A., Boynton, G.M., Fine, I.: A model of ganglion axon pathways accounts for percepts elicited by retinal implants. Scientific Reports 9(1), 1–16 (Jun 2019), https://doi.org/10.1038/s41598-019-45416-4
  * [5] Bratu, E., Dwyer, R., Noble, J.: A Graph-Based Method for Optimal Active Electrode Selection in Cochlear Implants. In: Martel, A.L., Abolmaesumi, P., Stoyanov, D., Mateus, D., Zuluaga, M.A., Zhou, S.K., Racoceanu, D., Joskowicz, L. (eds.) Medical Image Computing and Computer Assisted Intervention – MICCAI 2020. pp. 34–43. Lecture Notes in Computer Science, Springer International Publishing, Cham (2020), https://doi.org/10.1007/978-3-030-59716-0˙4
  * [6] Chakravorti, S., Noble, J.H., Gifford, R.H., Dawant, B.M., O’Connell, B.P., Wang, J., Labadie, R.F.: Further Evidence of the Relationship Between Cochlear Implant Electrode Positioning and Hearing Outcomes. Otology & Neurotology 40(5), 617–624 (Jun 2019), https://doi.org/10.1097/MAO.0000000000002204
  * [7] Erickson-Davis, C., Korzybska, H.: What do blind people “see” with retinal prostheses? Observations and qualitative reports of epiretinal implant users. PLOS ONE 16(2), e0229189 (Feb 2021), https://doi.org/10.1371/journal.pone.0229189, publisher: Public Library of Science
  * [8] Finley, C.C., Holden, T.A., Holden, L.K., Whiting, B.R., Chole, R.A., Neely, G.J., Hullar, T.E., Skinner, M.W.: Role of Electrode Placement as a Contributor to Variability in Cochlear Implant Outcomes. Otology & Neurotology 29(7), 920–928 (Oct 2008), https://doi.org/10.1097/MAO.0b013e318184f492
  * [9] Granley, J., Relic, L., Beyeler, M.: A Hybrid Neural Autoencoder for Sensory Neuroprostheses and Its Applications in Bionic Vision (May 2022), http://arxiv.org/abs/2205.13623, arXiv:2205.13623 [cs]
  * [10] Luo, Y.H.L., Zhong, J.J., Clemo, M., da Cruz, L.: Long-term Repeatability and Reproducibility of Phosphene Characteristics in Chronically Implanted Argus II Retinal Prosthesis Subjects. American Journal of Ophthalmology 170, 100–109 (Oct 2016), https://doi.org/10.1016/j.ajo.2016.07.021
  * [11] Nemhauser, G.L., Wolsey, L.A., Fisher, M.L.: An analysis of approximations for maximizing submodular set functions—I. Mathematical Programming 14(1), 265–294 (Dec 1978), https://doi.org/10.1007/BF01588971
  * [12] Rattay, F., Resatz, S.: Effective electrode configuration for selective stimulation with inner eye prostheses. IEEE Trans Biomed Eng 51(9), 1659–64 (Sep 2004), https://doi.org/10.1109/TBME.2004.828044
  * [13] Rizzo, J.F., Wyatt, J., Loewenstein, J., Kelly, S., Shire, D.: Perceptual Efficacy of Electrical Stimulation of Human Retina with a Microelectrode Array during Short-Term Surgical Trials. Investigative Ophthalmology & Visual Science 44(12), 5362–5369 (Dec 2003), https://doi.org/10.1167/iovs.02-0817
  * [14] Shah, N.P., Madugula, S., Grosberg, L., Mena, G., Tandon, P., Hottowy, P., Sher, A., Litke, A., Mitra, S., Chichilnisky, E.: Optimization of Electrical Stimulation for a High-Fidelity Artificial Retina. In: 2019 9th International IEEE/EMBS Conference on Neural Engineering (NER). pp. 714–718 (Mar 2019), https://doi.org/10.1109/NER.2019.8716987, iSSN: 1948-3546
  * [15] Wilke, R.G.H., Moghadam, G.K., Lovell, N.H., Suaning, G.J., Dokos, S.: Electric crosstalk impairs spatial resolution of multi-electrode arrays in retinal implants. Journal of Neural Engineering 8(4), 046016 (Jun 2011), https://doi.org/10.1088/1741-2560/8/4/046016