# Complex Locomotion Skill Learning
via Differentiable Physics

Yu Fang
University of Pennsylvania
<EMAIL_ADDRESS>
&Jiancheng Liu 11footnotemark: 1
<EMAIL_ADDRESS>
&Mingrui Zhang11footnotemark: 1
Imperial College London
<EMAIL_ADDRESS>
&Jiasheng Zhang
Taichi Graphics
<EMAIL_ADDRESS>
&Yidong Ma
Taichi Graphics
<EMAIL_ADDRESS>
&Minchen Li
University of California, Los Angeles
<EMAIL_ADDRESS>
&Yuanming Hu
Taichi Graphics
<EMAIL_ADDRESS>
&Chenfanfu Jiang
University of California, Los Angeles
<EMAIL_ADDRESS>
&Tiantian Liu
Taichi Graphics
<EMAIL_ADDRESS>
Equal contribution. The work was done when joint first authors were research
assistants at Taichi Graphics.Correspondence to Tiantian Liu

###### Abstract

Differentiable physics enables efficient gradient-based optimizations of
neural network (NN) controllers. However, existing work typically only
delivers NN controllers with limited capability and generalizability. We
present a practical learning framework that outputs unified NN controllers
capable of tasks with significantly improved complexity and diversity. To
systematically improve training robustness and efficiency, we investigated a
suite of improvements over the baseline approach, including periodic
activation functions, and tailored loss functions. In addition, we find our
adoption of batching and an Adam optimizer effective in training complex
locomotion tasks. We evaluate our framework on differentiable mass-spring and
material point method (MPM) simulations, with challenging locomotion tasks and
multiple robot designs. Experiments show that our learning framework, based on
differentiable physics, delivers better results than reinforcement learning
and converges much faster. We demonstrate that users can interactively control
soft robot locomotion and switch among multiple goals with specified velocity,
height, and direction instructions using a unified NN controller trained in
our system. Code is available at https://github.com/erizmr/Complex-locomotion-
skill-learning-via-differentiable-physics.

_Keywords_ Differentiable physics $\cdot$ locomotion skill learning $\cdot$
physical simulation

##  1 Introduction

Figure 1: Our learning system is robust and versatile, supporting various
simulation methods, robot designs, and locomotion tasks with continuously
controllable target velocity and heights.

Differentiable physical simulators deliver accurate analytical gradients of
physical simulations, opening up a promising stage for efficient neural
network (NN) controller training via gradient descent. Existing research
demonstrates that on simple tasks, learning systems via differentiable physics
can effectively leverage simulation gradient information and converge orders
of magnitude faster (see, e.g., de Avila Belbute-Peres et al. (2018); Hu et
al. (2019)) than reinforcement learning.

However, the capability of existing differentiable physics based learning
systems is relatively limited. Typically, optimized controllers can only
achieve relatively simple single-goal tasks (e.g., moving in one specific
direction, as in ChainQueen Hu et al. (2019)). Those learned controllers often
have difficulty generalizing the task to a perturbed version.

In this work, we propose a learning framework for complex locomotion skill
learning via differentiable physics. The complexity of our tasks comes from
two aspects: first, the degrees of freedom of our soft agents are
significantly higher compared with the rigid-body ones; second, the agents are
expected to learn multiple skills at the same time. We systematically propose
a suite of enhancements (Fig. 4) over existing training approaches (such as Hu
et al. (2020); Huang et al. (2021)) to significantly improve the efficiency,
robustness, and generalizability of learning systems based on differentiable
physics.

In addition, we investigated the contributions of each enhancements to
training in detail by a series of ablation studies. We also evaluated our
framework through comparisons against Proximal Policy Optimization (PPO)
Schulman et al. (2017), a state-of-the-art reinforcement learning algorithm.
Results show that our system is simple yet effective and has a much-improved
convergence rate compared to PPO.

We demonstrate the versatility of our framework via various physical
simulation environments (mass-spring systems and the moving least squares
material point method, MLS-MPM Hu et al. (2018)), robot designs (structured
and irregular), locomotion tasks (moving, jumping, turning around, all at a
controllable velocity). To the best of our knowledge, this is the first time
when a neural network can be trained via differentiable physics to achieve
tasks of such complexity (Fig. 1). Our agent is simultaneously trained with
multiple goals, the learned skills can be continuously interpolated. For
example, we can control both the velocity and height of an agent while it
runs. The trained agent can also be controlled in real-time, enabling direct
applications in soft robot control and video games. In addition, our method is
successfully applied to simple manipulation tasks (shown in Figure 13), which
indicate that our method has the potential to handle tasks beyond locomotion.
In summary, our key contributions are listed below:

  * •

To the best of our knowledge, we proposed the first differentiable physics
based learning framework for multiple locomotion skills learning on soft
agents using a single network.

  * •

We developed an end-to-end differentiable physical simulation environment for
deformable robot locomotion skills learning, which supports mass-spring system
and material point method as dynamic backends.

  * •

We systematically investigated the key factors contributing to the
differentiable physics based learning framework. We believe our investigation
can inspire further researches to explore the possibilities of differentiable
physics for more complex learning tasks.

##  2 Differentiable Simulation Environments

###  2.1 Simulation settings

Our learning framework supports multiple types of differentiable physically-
based deformable body simulators. The simulator can be seen as a black box
that outputs the states for the next time step given the actuation signal and
states from the current time step. The NN controller takes the information
given by the simulator as part of its input features and decides the actuation
signal for the simulator in the next time step. With our automatic
differentiation system, we can evaluate the gradient information to guide the
training of our NN controllers. We use the mass-spring systems and the moving
least squares material point method (MLS-MPM) as our simulation environments.
The design choices for each environment can be found in the appendix.

###  2.2 Agent structure and actuation mode

Figure 2: 3D agents collection. These agents are designed with simple stacked
cubes or complex handcrafted meshes.

Our framework supports a variety of differently shaped agents as shown in Fig.
2 and 3. Our actuation signal exerts forces along with the "muscle" directions
of the simulated agents. This signal falls within the interval
$\left[-1,+1\right]$, where its sign determines whether an agent wants to
contract (negative sign) or relax (positive sign) a muscle and its absolute
value determines the magnitude of the force. In mass-spring systems, the
muscle directions are represented by the spring directions, we change the
rest-length of springs to generate forces. The actuation signal is used to
scale the rest-length up to some limit (usually at $20\%$). In addition, only
activated springs, i.e. actuators (springs marked red or blue shown in Fig.
1), are able to generate forces according to control signals.

Figure 3: 2D agents collection.

In material point methods, the muscle directions are always along the vertical
axis in material space. In our implementation, we simply modify the Cauchy
stress in the vertical direction of material space to apply this force. Again,
the force is scaled to a user-defined bound which can be adjusted for
different simulations.

##  3 Learning Framework

Figure 4: Framework overview. Simulation instances are batched and executed in
parallel on GPUs. The whole system is end-to-end differentiable, and we use
gradient-based algorithms to optimize the neural network controller weights.
Each time step involves evaluating a NN controller inspired by SIREN Sitzmann
et al. (2020), and a differentiable simulator time integration implemented
using DiffTaichi Hu et al. (2020). Simulation states are fed back to the
initial state pool to improve the richness of training sets and thereby the
robustness of the resulted NN controller. A tailored loss function (as shown
in section 3.2) is designed for each task.

In this section, we describe our training framework based on differentiable
physics (Fig. 4). In summary, we develop a differentiable simulator with an
embedded NN controller. In each time step, the program performs an NN
controller inference and a differentiable simulator time integration. After a
few hundred time steps, the program back-propagates gradients end-to-end via
reverse-mode automatic differentiation, and updates the weights of the neural
network using the Adam optimizer Kingma and Ba (2014).

###  3.1 Task Representation

Given an agent $A$, its position and corresponding velocity are denoted as
${\bf{x}}=(x,y,z)$ and ${\bf{v}}=(u,v,w)$. We presents three tasks (running,
jumping and crawling) for 2D agents, two (running and rotating) for 3D. In a
whole simulation with total steps $\mathcal{T}$, an agent is expected to
achieve all goals in a goal sequence $\mathcal{G}=\\{{\bf G}_{1},{\bf
G}_{2}...{\bf G}_{n}|n\mathcal{P}=\mathcal{T}\\}$, where $\mathcal{P}$ is the
time period for an agent to achieve a goal, $n$ is the number of goals . For
each goal $\bf{G}$, there are multiple tasks in it. Each task is encoded using
a target value $g$, which instructs the controller to drive the agent. For
example, ${\bf G}=\\{g_{v},g_{h}\\}$ represent performing running and jumping
simultaneously, the target velocity and height are $g_{v}$ and $g_{h}$. The
agent is expected to switch between multiple goals during one simulation.

#### Running.

The running task is defined as a velocity tracking problem. Given a time step
$t$, the center of mass $\bf c$ is used to represent the position of the
agent, its velocity is defined as the averaged velocity $\tilde{\bf v}$ over a
time period $\mathcal{P}_{r}$, ($\mathcal{P}_{r}<\mathcal{P}$). Then the agent
is expected to achieve a velocity $\tilde{\bf v}$ close to the target velocity
$g_{v}$.

$\displaystyle\tilde{{\bf v}}(t)=\frac{1}{\mathcal{P}_{r}}\left({\bf
c}(t)-{\bf c}(t-\mathcal{P}_{r})\right)$ (1)

#### Jumping.

The jumping task is defined as reaching a given height $g_{h}$ in a time
period $\mathcal{P}$. Due to the gravity, it is not possible for a running
robot to stay away from the ground. Therefore, the jumping height
$\tilde{h}_{h}(t)$ of an agent is defined by the maximum vertical position of
its lowest point in a certain period $\mathcal{P}$. To be more formally, given
$\mathcal{{\bf S}}$ a set of nodes (or points) that constitute the agent, $h$
is a function that can extract the vertical position given a node (or point),
the jumping height over the $n$-th period is defined as

$\displaystyle\tilde{h}_{h}(n)=\max_{t\in[0,\mathcal{P}]}\min_{s\in\mathcal{{\bf
S}}}h(s,t+n\mathcal{P})$ (2)

where $n$ is the index of the period in a whole simulation.

#### Crawling.

The crawling task is defined as lowering the highest point as much as
possible. The target of crawling $g_{c}$ is defined as an indicator function.
The agent is controlled to crawl if $g_{c}=1$ otherwise if $g_{c}=0$. Since
crawling is a status that needs to be maintained, the crawling height is
defined as the vertical position of its highest point at each time step.

$\displaystyle\tilde{h}_{c}(t)=\max_{s\in\mathcal{{\bf S}}}h(s,t)$ (3)

###  3.2 Loss functions

Consider a locomotion task where an agent is instructed to move at a specified
velocity, we think a proper loss function should satisfy the following
requirements:

  1. 1.

Periodicity. The agent is expected to move periodically following a predefined
cyclic activation signal. Therefore, the loss function should encourage
periodic motions.

  2. 2.

Delayed evaluation. Due to inertia, it takes time for the agent to start
running or adjust running velocity. Therefore, an ideal loss function should
take this delay into consideration.

  3. 3.

Fluctuation tolerance. During a single running cycle, requiring the center of
mass to move at a constant velocity at each time step may induce highly
fluctuated losses. To avoid that, our loss function should be smooth during
one motion cycle.

#### Task loss.

For each task, we defined a tailored loss shown in equation 4 according to the
requirements above. The loss for running $\mathcal{L}_{v}$ is defined as the
accumulation of the difference between the target and the agent’s velocity
from start to current time step. For jumping, we define a sparse loss
$\mathcal{L}_{h}$, which only evaluate once for each period. The crawling loss
$\mathcal{L}_{c}$ is defined as the accumulation of the crawling height if the
loss is applicable.

$\displaystyle\mathcal{L}=\lambda_{v}\underbrace{\sum_{n\in\mathcal{T}}\sum_{t=\mathcal{P}_{r}}^{\mathcal{P}}(\tilde{{\bf
v}}(t)-g_{v}(n))^{2}}_{\mathcal{L}_{v}}+$
$\displaystyle\lambda_{h}\underbrace{\sum_{n\in\mathcal{T}}(\tilde{h}_{h}(n)-g_{h}(n))^{2}}_{\mathcal{L}_{h}}+$
$\displaystyle\lambda_{c}\underbrace{\sum_{n\in\mathcal{T}}\sum_{t=0}^{\mathcal{P}}g_{c}(n)\tilde{h}_{c}(t+n\mathcal{P})}_{\mathcal{L}_{c}}$
(4)

The $\lambda_{v},\lambda_{h},\lambda_{c}$ are weights for losses of dfferent
tasks. Our loss function accumulates the contributions of all steps in a
“sliding window”. This prevents the loss function from being too small as well
as the vanishing gradient problem. We find that having an accumulated loss
evaluated at every step works much better when we want to control the speed of
the agent explicitly.

#### Regularization.

Since the agents may keep shaking when no target velocity or height are given,
we introduce an actuation loss as a regularization term. The intuition is to
penalize comparatively large actuations when small goals value is given:

$\displaystyle\mathcal{L}_{a}=\sum_{n\in\mathcal{T}}\sum_{t\in[0,\mathcal{P}]}({\mathcal{A}}(t)-\mu||g_{v}(n)||)^{2}$
(5)

where $\mathcal{A}(t)$ is the actuation output of NN, $\mu$ is a normalizer.

###  3.3 Network architecture

We use two fully connected (FC) layers with sine activation function as the
neural network. The network input vector consists of three parts:

  1. 1.

Periodic control signal of the same period with different phases, to encourage
periodic actuation. In this work, we use sin waves as the periodic control
signal.

  2. 2.

State feature vector extracted from the current state of the agent. Using the
positions and velocities of all the vertices of the simulated object as the
input feature is not a good idea since these quantities are neither
translation-invariant nor rotation-invariant. We use a “centralized pose” to
remove global translation and rotation information. To be more specific, we
subtract the position of the center of mass (CoM) in both 2D and 3D cases to
remove the global translation. In 3D cases, we also compute the rotation
around the vertical axis of the agent to remove its global rotation.

  3. 3.

Targets that encode the task information such as running velocity, jumping
height, and orientation for 3D agents. Different tasks use different channels
to represents their target values. During the training, we assign random
values to these targets. When validating the model, we always test the agent
on a fixed set of target values. In interactive settings, we assign various
target values to control the agent’s motion. To amplify the importance of the
targets, we duplicate these channels multiple times, try to make them
comparably important among the other input features.

###  3.4 Training

In the training process, we run one simulation in one training iteration. The
duration of the simulation is divided into several periods $\mathcal{P}$ as
mentioned in the section 3.1. For each $\mathcal{P}$, a goal $\bf G$ contains
several tasks is assigned. For example, ${\bf G}=\\{g_{v},g_{h}\\}$ represents
that the agent is expected to perform both running and jumping in this period.
The goal and its target values are generated randomly in a uniform
distribution. They are kept fixed inside one period $\mathcal{P}$, but varied
between different periods.

##  4 Experiments and Analysis

In this section, we systematically investigate the key factors that contribute
to the training, evaluate our framework through ablation studies and
comparisons against PPO Schulman et al. (2017), a state-of-the-art
reinforcement learning algorithm. All experiments in the ablation studies are
performed at least three times. The results of ablation studies are summarized
in Table 1 and Fig. 5.

Name | Training Setting | Norm. Valid. Loss (Alpaca)
---|---|---
OP | AF | BS | PS | SV | TG | LD | Run. | Jump. | Task
Full | Adam | sin | 32 | ✓ | ✓ | ✓ | ✓ | 0.58$\pm$0.01 | 0.79$\pm$0.01 | 0.74$\pm$0.01
Difftaichi* | SGD | tanh | 1 | ✓ | ✓ | ✓ | ✓ | 0.98$\pm$0.00 | 0.95$\pm$0.01 | 0.96$\pm$0.01
Full-BS | Adam | sin | 1 | ✓ | ✓ | ✓ | ✓ | 0.82$\pm$0.02 | 0.85$\pm$0.01 | 0.84$\pm$0.02
Full-OP | SGD | sin | 32 | ✓ | ✓ | ✓ | ✓ | 0.97$\pm$0.01 | 0.98$\pm$0.01 | 0.98$\pm$0.01
Full-AF | Adam | tanh | 32 | ✓ | ✓ | ✓ | ✓ | 0.67$\pm$0.02 | 0.81$\pm$0.01 | 0.77$\pm$0.02
Full-PS | Adam | sin | 32 | ✗ | ✓ | ✓ | ✓ | 0.97$\pm$0.00 | 0.99$\pm$0.00 | 0.99$\pm$0.00
Full-SV | Adam | sin | 32 | ✓ | ✗ | ✓ | ✓ | 0.54$\pm$0.01 | 0.84$\pm$0.01 | 0.77$\pm$0.01
Full-TG | Adam | sin | 32 | ✓ | ✓ | ✗ | ✓ | 0.99$\pm$0.00 | 0.95$\pm$0.00 | 0.97$\pm$0.00
Full-LD | Adam | sin | 32 | ✓ | ✓ | ✓ | ✗ | 0.99$\pm$0.01 | 0.74$\pm$0.01 | 0.82$\pm$0.01
Table 1:  Ablation summary on agent _Alpaca_. In this table, we show the
_Task_ , _Run_ and _Jump_ validation loss of proposed method and its ablated
versions. _Full_ represents the proposed method. The abbreviations under
training setting represent name of components to remove. To be more specific,
OP: Optimizer, AF: Activation Function, BS: Batch Size, PS: Periodic Signal,
SV: State Vector, TV: Targets, LD: Loss Design. *Difftaichi is an
implementation of Hu et al. (2020), enhanced with our loss design. It can be
observed that the _Full_ method achieves the best performance on _Task_ loss
among all models.

Figure 5: Summary of the ablation study. Here we show the summary of ablation
study on agent _Alpaca_. Each iteration indicates one training iteration,
which is composed of 1000 steps of physical simulation and one step of update
on weights of the controller. The result show that the _Full_ method achieves
the best performance. For more results of ablation studies on other agents,
please check the appendix.

###  4.1 Validation Metric

In validation, we test all the trained agents on a fixed set of goals. These
goals are represented by combination of targets targets values, which are
uniformly sampled by interpolating between the lower and upper bound of all
targets appeared in training.

###  4.2 Differentiable Physics Gradient Analysis

One concern in differentiable physics based learning is the that the gradient
may explode during the back-propagation through long simulation steps. Here we
visualize the distribution of gradients norm of a whole training process in
Fig. 6. It can be observed that most values of the gradients are distributed
inside the interval [-10, 10] in log scale, which indicates that our learning
approach provides stable gradients during training.

Figure 6: Gradient Analysis. The plots show the gradient distribution of
different agent. The x-axis represent the sum of gradients norm. The values
are drawn in log scale for a clear visualization. We record the sum of
gradients norm of 10000 iterations for one training. For each agent, the
experiments are performed at least three times, i.e., there are at least 30000
gradient samples for one agent.

#### Adam vs SGD.

Hu et al. (2020) shows that stochastic gradient descent (SGD) with
differentiable physics can achieve satisfying results for tasks with single
goals, e.g. running along one direction. However, this no longer holds for
tasks with multiple goals. To investigate the differences, we perform a
learning task with multiple goals on a 2D mass spring system. In this
experiment shown in the left of Fig. 5, SGD (curve Full-OP) makes minor
progress while a Adam optimizer with reduced momentum drops the validation
obviously within 10000 iterations. One key difference between Adam and SGD is
that Adam utilises the momentum while vanilla SGD doesn’t. To further
investigate the the effectiveness of the momentum hyperparameters $\beta$, we
perform a grid search for both $\beta_{1}$ and $\beta_{2}$ whose default
values are $0.9$ and $0.999$ respectively. The results show that reducing
$\beta_{2}$ to $0.9$ delivers us the best performance. We therefore use this
Adam optimizer with reduced momentum in all our experiments.

###  4.3 Learning techniques

#### Batching.

We adopt batching into differentiable physics and perform a series of
experiments as shown in the right of Fig. 5 to verify its effectiveness, with
batch size 1, 8, 32, and 64. It can be observed that without batching, there
is a large variance in training loss and it hardly converges. Batching with a
proper size helps reduce the variance and effectively improve the training
performance. We also find that an overlarge batching size does not help
training obviously. For example, batching with size 64 only improve the
performance marginally.

#### Loss design.

We compare our tailored loss function with a naive loss function. The naive
loss function defines the velocity as the position difference of the center of
mass between consecutive time steps, i.e., requiring the agent to move at a
constant velocity at each time step. The experiments show that the agent under
the naive loss design can not be properly trained to run, i.e., there are
almost no progress in running loss dropping.

Figure 7: Activation functions ablation study. The figure show the validation
loss between choices of different activation functions. From left to right,
the subplot shows the results of _Task_ , _Run_ and _Jump_ loss, respectively.

#### Activation function.

To validate the effectiveness of our $sin$ activation function, we replaced it
for hidden layer with various popular activation functions $tanh$, $relu$,
$gelu$ and $sigmoid$. For the output layer, the output actuation is expected
to be limited in the range -1 to 1 due to the physical constraints of an
actuator, e.g., a spring should not be overly compressed or stretched.
Therefore we made some modifications for the activation functions in the
output layer. To be more specific, we clamp the values larger than 1 or
smaller than -1. In addition, for activation functions whose values are
constantly zero (or very small) when input is negative, we modify the $relu$
and $gelu$ by subtracting the value by 1, for $sigmoid$ we map the its value
from [0, 1] to [-1, 1]. The results in Fig.7 show that using $sin$ as
activation functions gives the best performance among all three losses. For
more results of activation functions ablation studies, please check the table
6 in appendix.

###  4.4 Importance of input features

Figure 8: Input features ablation study. Full-PS, Full-SV and Full-TG indicate
the models trained without Periodic Signal, State Vector or Targets
respectively.

#### Ablation on input features.

The input features consist of three parts: periodic control signal, state
vector, and targets, as mentioned in section 3.3. We investigate each part of
the contribution to training by ablation studies as shown in Fig.8. The
results show that each of the three parts serves as a necessary component for
a converged training. The periodic control signal serves as an important role.
The agents can hardly move without this signal as the corresponding task loss
does not drop along the training (shown by the red curves in Fig.8). The
targets are another essential part, without which the optimizer can hardly
work after a minor progress as shown by the yellow curves in Fig.8. Although
the state vector is less important in terms of reducing the loss, the overall
task performance is degraded without the state vector as shown by the green
curves in Fig.8.

###  4.5 Benchmark against reinforcement learning

Figure 9: Comparison on different agents. Both our method and PPO run on GPUs.
The solid and dashed lines show the validation loss of our method and PPO,
respectively. The agent design is also shown in the figure. Springs marked in
red are actuators.

We implement the standard proximal policy optimization (PPO) benchmark in
multiple goals settings as shown in Fig. 9. The PPO agent can make progress in
single goal learning and learn certain locomotion patterns to move toward the
goal. Yet as the task gets complicated, PPO often gets trapped in the local
minima. We try our best to fine-tune the PPO hyper-parameters, but still find
it tends to overfit to one goal but fails to find a balance between different
goals. It can be observed in Fig. 9 where the agents trained by PPO can learn
to jump but struggle to run.

##  5 Related work

#### Differentiable simulation.

Differentiable physical simulation is getting increasingly more attention in
the learning community. Two families of methods exist: the first family uses
neural networks to approximate physical simulators (Battaglia et al., 2016;
Chang et al., 2016; Mrowca et al., 2018; Li et al., 2018). Differentiating the
approximating NNs then yields gradients of the (approximate) physical
simulation. The second family is more accurate and direct: many of these
methods using differentiable programming (specifically, reverse-mode automatic
differentiation) to implement physical simulators (Degrave et al., 2016; de
Avila Belbute-Peres et al., 2018; Schenck and Fox, 2018; Heiden et al., 2019;
Hu et al., 2019, 2020; Huang et al., 2021). Automatic differentiation works
well for explicit time integrators, but when it comes to implicit time
integration, people often adopt the adjoint methods (Bern et al., 2019;
Geilinger et al., 2020), LCP (de Avila Belbute-Peres et al., 2018) and QR
decompositions (Liang et al., 2019; Qiao et al., 2020). In this work, we
leverage DiffTaichi (Hu et al., 2020), an automatic differentiation system, to
create high-performance parallel differentiable simulators and the built-in NN
controllers.

#### Locomotion skill learning.

Producing physically plausible characters with various locomotion skills is a
challenging problem. One classic approach is to manually design locomotion
controllers subject to physics laws of motion in order to generate walking Ye
and Liu (2010); Coros et al. (2010) or bicycling Tan et al. (2014) characters.
These physically-based controllers often rely on a complex set of control
parameters and are difficult to design and time-consuming to optimize. Another
line of research suggests using data to control characters kinematically.
Given a set of data clips, controllers can be made to select the best fit clip
to properly react to certain situations Safonova and Hodgins (2007); Lee et
al. (2010); Liu et al. (2010). These kinematic models use the motion clips to
build a state machine and add transitions between similar frames in adjacent
states. Although being able to produce higher quality motions than most
simulation-based methods, the kinematic methods lack the ability to synthesize
behaviors for unseen situations. With the recent advances in deep learning
techniques, attempts have been explored to incorporate reinforcement learning
(RL) into locomotion skill learning Peng et al. (2015); Liu and Hodgins
(2017); Peng et al. (2018); Park et al. (2019). These modern controllers gain
their effectiveness either from tracking high-quality reference motion clips
or from cleverly designed rewards to imitate the reference. However, the
exploration space for RL is usually prohibitively large to achieve complicated
target motions. Carefully designed early termination strategies Ma et al.
(2021); Won et al. (2020), better optimization methods Yang and Yin (2021),
and adversarial RL schemes Peng et al. (2021) enable these RL-based methods to
achieve richer behaviors. However, it is still time-consuming to get a well-
trained RL model in complex tasks such as soft robot control. Our method, on
the other hand, leverages the differentiable simulation framework and can
achieve much better convergence behavior compared to RL-based methods.

##  6 Limitations and Conclusions

Currently, difficulty of the learning tasks depends on robot design and
physical parameters. The training performance may be degraded with improper
robot designs. For instance, in our stool case where unnecessary actuators on
its body are allowed, it is more difficult to achieve good results. Offloading
the physical parameter tuning and robot designing to an automatic pipeline
will be an interesting future research direction.

To summarize, we have presented an effective end-to-end differentiable physics
based learning framework for soft robot complex locomotion skills learning. We
systematically enhance classical differentiable physics learning systems with
a suite of techniques and investigated the key factors contributing to the
training in detail. We show that our framework can provide stable gradients
without explosion during a simulation with hundreds of steps. Together with
batching, the Adam with reduced momentum dramatically outperforms the
stochastic gradient descent (SGD) on complex tasks with multiple goals.
Benefited by our tailored loss functions, the network take advantages of the
three parts of the input features. We found that the periodic control signal
dominates the actuation signal, while state vector and task goals have weaker
effects. To penalize the comparatively large actuations when small goals value
is given, we additionally pose an actuation loss as a regularization term. It
can effectively suppress the noises induced by the overly large inputs, which
help the agents move more smoothly. We also compare our method with the state-
of-the-art reinforcement learning proximal policy optimization. Our method
shows advantages in both training performance and convergence efficiency on
tasks with multiple goals. In addition, we successfully applied the proposed
method to simple manipulation tasks, which indicates the method has the
potential to generalize to tasks beyond locomotion.

Our framework enables users to flexibly design robots and to teach them
locomotion skills. The trained agents can be manipulated to smoothly switch
locomotion tasks such as running, jumping, and crawling with different speeds
and orientations, interactively. To the best of our knowledge, this is the
first time a learning system based on differentiable physics can deliver
controllers with such robustness, flexibility, and efficiency.

## References

  * de Avila Belbute-Peres et al. [2018] Filipe de Avila Belbute-Peres, Kevin Smith, Kelsey Allen, Josh Tenenbaum, and J Zico Kolter.  End-to-end differentiable physics for learning and control.  In _Advances in Neural Information Processing Systems_ , pages 7178–7189, 2018.
  * Hu et al. [2019] Yuanming Hu, Jiancheng Liu, Andrew Spielberg, Joshua B Tenenbaum, William T Freeman, Jiajun Wu, Daniela Rus, and Wojciech Matusik.  Chainqueen: A real-time differentiable physical simulator for soft robotics.  In _2019 International Conference on Robotics and Automation (ICRA)_ , pages 6265–6271. IEEE, 2019.
  * Hu et al. [2020] Yuanming Hu, Luke Anderson, Tzu-Mao Li, Qi Sun, Nathan Carr, Jonathan Ragan-Kelley, and Frédo Durand.  Difftaichi: Differentiable programming for physical simulation.  _ICLR_ , 2020.
  * Huang et al. [2021] Zhiao Huang, Yuanming Hu, Tao Du, Siyuan Zhou, Hao Su, Joshua B Tenenbaum, and Chuang Gan.  Plasticinelab: A soft-body manipulation benchmark with differentiable physics.  _arXiv preprint arXiv:2104.03311_ , 2021.
  * Schulman et al. [2017] John Schulman, Filip Wolski, Prafulla Dhariwal, Alec Radford, and Oleg Klimov.  Proximal policy optimization algorithms.  _arXiv preprint arXiv:1707.06347_ , 2017.
  * Hu et al. [2018] Yuanming Hu, Yu Fang, Ziheng Ge, Ziyin Qu, Yixin Zhu, Andre Pradhana, and Chenfanfu Jiang.  A moving least squares material point method with displacement discontinuity and two-way rigid body coupling.  _ACM Transactions on Graphics (TOG)_ , 37(4):1–14, 2018.
  * Sitzmann et al. [2020] Vincent Sitzmann, Julien Martel, Alexander Bergman, David Lindell, and Gordon Wetzstein.  Implicit neural representations with periodic activation functions.  _Advances in Neural Information Processing Systems_ , 33, 2020.
  * Kingma and Ba [2014] Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ , 2014.
  * Battaglia et al. [2016] Peter Battaglia, Razvan Pascanu, Matthew Lai, Danilo Jimenez Rezende, et al.  Interaction networks for learning about objects, relations and physics.  In _Advances in neural information processing systems_ , pages 4502–4510, 2016.
  * Chang et al. [2016] Michael B Chang, Tomer Ullman, Antonio Torralba, and Joshua B Tenenbaum.  A compositional object-based approach to learning physical dynamics.  _ICLR_ , 2016.
  * Mrowca et al. [2018] Damian Mrowca, Chengxu Zhuang, Elias Wang, Nick Haber, Li Fei-Fei, Joshua B Tenenbaum, and Daniel LK Yamins.  Flexible neural representation for physics prediction.  _arxiv preprint arXiv:1806.08047_ , 2018.
  * Li et al. [2018] Yunzhu Li, Jiajun Wu, Russ Tedrake, Joshua B Tenenbaum, and Antonio Torralba.  Learning particle dynamics for manipulating rigid bodies, deformable objects, and fluids.  _arXiv preprint arXiv:1810.01566_ , 2018.
  * Degrave et al. [2016] Jonas Degrave, Michiel Hermans, Joni Dambre, et al.  A differentiable physics engine for deep learning in robotics.  _arXiv preprint arXiv:1611.01652_ , 2016.
  * Schenck and Fox [2018] Connor Schenck and Dieter Fox.  Spnets: Differentiable fluid dynamics for deep neural networks.  _arXiv preprint arXiv:1806.06094_ , 2018.
  * Heiden et al. [2019] Eric Heiden, David Millard, Hejia Zhang, and Gaurav S Sukhatme.  Interactive differentiable simulation.  _arXiv preprint arXiv:1905.10706_ , 2019.
  * Bern et al. [2019] James M Bern, Pol Banzet, Roi Poranne, and Stelian Coros.  Trajectory optimization for cable-driven soft robot locomotion.  In _Robotics: Science and Systems_ , 2019.
  * Geilinger et al. [2020] Moritz Geilinger, David Hahn, Jonas Zehnder, Moritz Bächer, Bernhard Thomaszewski, and Stelian Coros.  Add: Analytically differentiable dynamics for multi-body systems with frictional contact.  _arXiv preprint arXiv:2007.00987_ , 2020.
  * Liang et al. [2019] Junbang Liang, Ming C Lin, and Vladlen Koltun.  Differentiable cloth simulation for inverse problems.  _Advances in Neural Information Processing Systems_ , 2019.
  * Qiao et al. [2020] Yi-Ling Qiao, Junbang Liang, Vladlen Koltun, and Ming C Lin.  Scalable differentiable physics for learning and control.  _arXiv preprint arXiv:2007.02168_ , 2020.
  * Ye and Liu [2010] Yuting Ye and C Karen Liu.  Optimal feedback control for character animation using an abstract model.  In _ACM SIGGRAPH 2010 papers_ , pages 1–9. 2010.
  * Coros et al. [2010] Stelian Coros, Philippe Beaudoin, and Michiel Van de Panne.  Generalized biped walking control.  _ACM Transactions On Graphics (TOG)_ , 29(4):1–9, 2010.
  * Tan et al. [2014] Jie Tan, Yuting Gu, C Karen Liu, and Greg Turk.  Learning bicycle stunts.  _ACM Transactions on Graphics (TOG)_ , 33(4):1–12, 2014.
  * Safonova and Hodgins [2007] Alla Safonova and Jessica K Hodgins.  Construction and optimal search of interpolated motion graphs.  In _ACM SIGGRAPH 2007 papers_ , pages 106–es. 2007.
  * Lee et al. [2010] Yongjoon Lee, Kevin Wampler, Gilbert Bernstein, Jovan Popović, and Zoran Popović.  Motion fields for interactive character locomotion.  In _ACM SIGGRAPH Asia 2010 papers_ , pages 1–8. 2010.
  * Liu et al. [2010] Libin Liu, KangKang Yin, Michiel van de Panne, Tianjia Shao, and Weiwei Xu.  Sampling-based contact-rich motion control.  In _ACM SIGGRAPH 2010 papers_ , pages 1–10. 2010.
  * Peng et al. [2015] Xue Bin Peng, Glen Berseth, and Michiel Van de Panne.  Dynamic terrain traversal skills using reinforcement learning.  _ACM Transactions on Graphics (TOG)_ , 34(4):1–11, 2015.
  * Liu and Hodgins [2017] Libin Liu and Jessica Hodgins.  Learning to schedule control fragments for physics-based characters using deep q-learning.  _ACM Transactions on Graphics (TOG)_ , 36(3):1–14, 2017.
  * Peng et al. [2018] Xue Bin Peng, Pieter Abbeel, Sergey Levine, and Michiel van de Panne.  Deepmimic: Example-guided deep reinforcement learning of physics-based character skills.  _ACM Transactions on Graphics (TOG)_ , 37(4):1–14, 2018.
  * Park et al. [2019] Soohwan Park, Hoseok Ryu, Seyoung Lee, Sunmin Lee, and Jehee Lee.  Learning predict-and-simulate policies from unorganized human motion data.  _ACM Transactions on Graphics (TOG)_ , 38(6):1–11, 2019.
  * Ma et al. [2021] Li-Ke Ma, Zeshi Yang, Xin Tong, Baining Guo, and KangKang Yin.  Learning and exploring motor skills with spacetime bounds.  _arXiv preprint arXiv:2103.16807_ , 2021.
  * Won et al. [2020] Jungdam Won, Deepak Gopinath, and Jessica Hodgins.  A scalable approach to control diverse behaviors for physically simulated characters.  _ACM Transactions on Graphics (TOG)_ , 39(4):33–1, 2020.
  * Yang and Yin [2021] Zeshi Yang and Zhiqi Yin.  Efficient hyperparameter optimization for physics-based character animation.  _Proceedings of the ACM on Computer Graphics and Interactive Techniques_ , 4(1):1–19, 2021.
  * Peng et al. [2021] Xue Bin Peng, Ze Ma, Pieter Abbeel, Sergey Levine, and Angjoo Kanazawa.  Amp: Adversarial motion priors for stylized physics-based character control.  _arXiv preprint arXiv:2104.02180_ , 2021.
  * Kostrikov [2018] Ilya Kostrikov.  Pytorch implementations of reinforcement learning algorithms.  https://github.com/ikostrikov/pytorch-a2c-ppo-acktr-gail, 2018.
  * Baraff and Witkin [1998] David Baraff and Andrew Witkin.  Large steps in cloth simulation.  In _Proceedings of the 25th annual conference on Computer graphics and interactive techniques_ , pages 43–54, 1998.
  * Jiang et al. [2015] Chenfanfu Jiang, Craig Schroeder, Andrew Selle, Joseph Teran, and Alexey Stomakhin.  The affine particle-in-cell method.  _ACM Transactions on Graphics (TOG)_ , 34(4):1–10, 2015.

##  Appendix A Reinforcement learning Setting

#### Environment.

We use the open-source implementation of PPO (Kostrikov [2018]) in our
environments. Part of important hyper-parameters are listed in the table 2.

Table 2: PPO hyper-parameters Parameter | Values
---|---
learning rate | $2.5e-4$
entropy coef | 0.01
value loss coef | 0.5
number of processes | 8
number of simulation steps. | 1000

#### Reward Design.

We design a reward function for PPO based on our loss functions. Recall the
velocity part of equation 4. We define the total velocity reward as:

$\displaystyle\mathcal{R}_{v}$
$\displaystyle=\lambda_{v}\sum_{n\in\mathcal{T}}\sum_{t\in[\mathcal{P}_{r},\mathcal{P}]}g_{v}(n)^{2}-(\tilde{{\bf
v}}(t)-g_{v}(n))^{2}$ (6)

If we directly split the reward into each time steps by $t$, the rewards of
first $\mathcal{P}_{r}$ time steps would be zero and it is too difficult for
PPO to learn the policy. Therefore, we modified the reward function to tackle
this issue. By equation 1, we have

$\displaystyle\mathcal{R}_{v}$
$\displaystyle=\lambda_{v}\sum_{n\in\mathcal{T}}\sum_{t\in[\mathcal{P}_{r},\mathcal{P}]}g_{v}(n)^{2}-\left[\frac{1}{\mathcal{P}_{r}}\left({\bf
c}(t)-{\bf c}(t-\mathcal{P}_{r})\right)-g_{v}(n)\right]^{2}$ (7)
$\displaystyle=\frac{\lambda_{v}}{\mathcal{P}_{r}^{2}}\sum_{n\in\mathcal{T}}\sum_{t\in[\mathcal{P}_{r},\mathcal{P}]}[\mathcal{P}_{r}g_{v}(n)]^{2}-\left[{\bf
c}(t)-({\bf c}(t-\mathcal{P}_{r})+\mathcal{P}_{r}g_{v}(n))\right]^{2}.$ (8)

We can define a function
$f_{n}(t^{\prime},t)=[\mathcal{P}_{r}g_{v}(n)]^{2}-\left[{\bf
c}(t^{\prime})-({\bf
c}(t-\mathcal{P}_{r})+\mathcal{P}_{r}g_{v}(n))\right]^{2}$ and substitute it
back to equation 8. The reward can be re-written as

$\displaystyle\mathcal{R}_{v}$
$\displaystyle=\frac{\lambda_{v}}{\mathcal{P}_{r}^{2}}\sum_{n\in\mathcal{T}}\sum_{t\in[\mathcal{P}_{r},\mathcal{P}]}f_{n}(t,t)$
(9)

Refer to that for any $t\in\mathcal{P}$, we have
$f_{n}(t-\mathcal{P}_{r},t)=0$, which means

$\displaystyle\mathcal{R}_{v}$
$\displaystyle=\frac{\lambda_{v}}{\mathcal{P}_{r}^{2}}\sum_{n\in\mathcal{T}}\sum_{t\in[\mathcal{P}_{r},\mathcal{P}]}\sum_{\Delta
t\in[0,\mathcal{P}_{r})}f_{n}(t-\Delta t,t)-f_{n}(t-\Delta t-1,t)$ (10)

Let $t^{\prime}=t-\Delta t$

$\displaystyle\mathcal{R}_{v}=$
$\displaystyle\frac{\lambda_{v}}{\mathcal{P}_{r}^{2}}\sum_{n\in\mathcal{T}}\sum_{t^{\prime}\in[0,\mathcal{P}]}\sum_{\Delta
t=0}^{\min\left(t^{\prime},\mathcal{P}_{r}\right)}f_{n}(t^{\prime},t^{\prime}+\Delta
t)-f_{n}(t^{\prime}-1,t^{\prime}+\Delta t)$ (11) $\displaystyle=$
$\displaystyle\frac{\lambda_{v}}{\mathcal{P}_{r}^{2}}\sum_{n\in\mathcal{T}}\sum_{t^{\prime}\in[0,\mathcal{P}]}\sum_{\Delta
t=0}^{\min\left(t^{\prime},\mathcal{P}_{r}\right)}-\left[{\bf
c}(t^{\prime})-({\bf c}(t^{\prime}+\Delta
t-\mathcal{P}_{r})+\mathcal{P}_{r}g_{v}(n))\right]^{2}+$
$\displaystyle\left[{\bf c}(t^{\prime}-1)-({\bf c}(t^{\prime}+\Delta
t-\mathcal{P}_{r})+\mathcal{P}_{r}g_{v}(n))\right]^{2}$ (12)

So we can split the whole reward into each time step $t$ as:

$\displaystyle\mathcal{R}_{v}(t)=$ $\displaystyle\sum_{\Delta
t=0}^{\min\left(t,\mathcal{P}_{r}\right)}\left[{\bf c}(t-1)-({\bf c}(t+\Delta
t-\mathcal{P}_{r})+\mathcal{P}_{r}g_{v}(n))\right]^{2}-$
$\displaystyle\left[{\bf c}(t)-({\bf c}(t+\Delta
t-\mathcal{P}_{r})+\mathcal{P}_{r}g_{v}(n))\right]^{2}$ (13)

Recall that the goal of a jumping task is to reach a specified height goal
$g_{h}$, the intuition for jumping reward is to encourage the agent to improve
its maximum height until reaching the height.

$\displaystyle\mathcal{R}_{h}(t)=\lambda_{h}\sum_{n\in\mathcal{T}}(\min_{s\in\mathcal{{\bf
S}}}h(s,t+n\mathcal{P})-g_{h}(n))^{2}-(\tilde{h}_{h}(n)-g_{h}(n))^{2}$ (14)

where the $\min_{s\in\mathcal{{\bf S}}}h(s,t+n\mathcal{P})$ is the agent
height at current time step as mentioned in 3.1 and $\tilde{h}_{h}(n)$ is the
max height record during one task period.

##  Appendix B Network Architecture

Figure 10: Network architecture.

Our network architecture is shown in Fig. 10. It is a two layers fully
connected network with a 64 channels hidden layer. The number of input and
output channels varies according to different agent designs.

##  Appendix C Simulation Setting

#### Mass-spring systems

We adopt the classic Hookean spring model to represent the elastic force and
use dashpot damping Baraff and Witkin [1998] as the damping force to simulate
mass-spring system. We found that the drag damping model used by Hu et al. Hu
et al. [2020] damps the gradient of the NN controller while the dashpot
damping model is able to generate vividly changing gradients. Adopting dashpot
damping model makes our agent more flexible.

#### Material point methods

We use MLS-MPM Hu et al. [2018] as our material point method simulator. We
further applied the affine particle-in-cell method Jiang et al. [2015] to
reduce the artificial damping. Due to the nature of MPM as a hybrid
Lagrangian-Eulerian method, it always requires a background grid during the
simulation. Naive MPM implementations fix the position of the background grid,
hence limit the moving range of the agent. We apply a dynamic background grid
that follows the agent all the time to overcome this problem.

##  Appendix D Interactive control

Figure 11: Left: Our system allows the user to control the agent’s motion
interactively. Right: the trajectory of a walking quadruped agent under user
control.

Smooth interpolation between different tasks is a key advantage of our
approach. In interactive settings, users can smoothly control the agent’s
motion via input devices such as a keyboard as shown in Fig. 11. Please refer
to our supplemental video for more details.

##  Appendix E 3D results

In addition to 2D cases, our framework can be applied to 3D cases seamlessly.
In a 3D space, an agent is additionally expected to control its orientation on
a plane, i.e., rotate. We design several 3D agents (shown in Fig. 2) and train
them to run and rotate. These trained agents can achieve specified goals given
control signals after trained for several hundred to a few thousand
iterations. Please check the supplemental videos for more details.

#### Loss function in 3D

The agents in 3D are trained with the loss function below. The loss for
running $\mathcal{L}_{v}$ is inherited from 2D settings. The rotation loss
$\mathcal{L}_{r}$ is defined as the accumulation of centralized point-wise
distance between current position and target rotated position.

$\displaystyle\mathcal{L}=\lambda_{v}\underbrace{\sum_{n\in\mathcal{T}}\sum_{t=\mathcal{P}_{r}}^{\mathcal{P}}(\tilde{{\bf
v}}(t)-g_{v}(n))^{2}}_{\mathcal{L}_{r}}+\lambda_{r}\underbrace{\sum_{n\in\mathcal{T}}\sum_{t=\mathcal{P}_{r}}^{\mathcal{P}}\sum_{s\in{\bf
S}}(({\bf s}(t)-{\bf c}(t))-{\bf R}({\bf s}(t-\mathcal{P}_{r})-{\bf
c}(t-\mathcal{P}_{r})))^{2}}_{\mathcal{L}_{r}}$ (15) $\displaystyle\bf
R=\begin{bmatrix}\cos(\mathcal{P}_{r}g_{\omega})&0&-\sin(\mathcal{P}_{r}g_{\omega})\\\
0&1&0\\\ \sin(\mathcal{P}_{r}g_{\omega})&0&\cos(\mathcal{P}_{r}g_{\omega})\\\
\end{bmatrix},$ (16)

where $g_{\omega}$ is the target angular velocity assembled in input features.
The format of $\bf R$ is derived by rotation matrix along XZ-plane.

Note, for solving running tasks, the rotation loss also plays an important
role. A significant issue of 3D agents running is accumulated orientation
error. Utilizing identity rotation matrix $\bf R$ helps adjust agent’s
posture. In turn’s of rotation tasks, the rotation center should be maintained
still which is equivalent to zero velocity. These two losses complement with
each other to achieve the final goals.

#### Friction model.

In real world, agents running is driven by static friction between "feet" and
ground. We configure and train agents with different contact models in our
experiments. With zero friction applied, training loss shows no evidence
decreasing and the agent is not able to move no matter what actions patterns
it learns. With a fully sticky surface, where only upwards vertical velocity
is kept and all other components of velocities are projected to zero after
contact, the 2D agents work well. However, in 3D, it turns out that the agent
tends to stuck. In practice, we adopt the classic slip with friction model
from physically-based simulation which involves both kinetic friction and
static friction. These two categories of contact are determined by how much
pressure force is applied on the surface and can be treated like a branched
function numerically which can be handled by our auto-differential system.
Please see the supplemental video for a visual comparison between slip and
sticky friction models.

Figure 12: Left: Training loss for zero friction, .4 friction and sticky
contact models. Middle: Agent suffers from sticky surface and cannot move
further. Right: Agent moves smoothly towards right.

##  Appendix F Generalization to Manipulation Tasks

Here we provide results of applying our method on simple manipulation tasks.
In these tasks, the agent is expected to manipulate an object (marked in
purple) to hit a target point (marked in green). This task has similar
periodicity property like locomotion. We designed two scenarios: _Juggle_ and
_Dribbble and Shot_ shown in Fig.13. The visual results are also shown in the
video in supplemental materials.

#### Juggle

In this scenario, the target point appears in the sky above the the agent. The
agent is expected to ’juggle’ the object to the target point.

#### Dribble and Shot

In this scenario, the target point keeps moving toward right. The agent is
expected to carry the object and shot it to the target point.

Figure 13: Manipulation tasks showcase. The upper plot shows the snapshots of
scenario ’Juggle’ and the lower one shows that of the ’Dribble and Shot’.
Figure 14: Training curve of the ’Juggle’ task. The loss is defined as the
distance between center of the object and the target point.

The experiments indicate that our method has the possibility to generalize to
tasks beyond locomotion.

##  Appendix G Detailed Results on Further Experiments

Here we show more detailed results on further experiments, including ablation
studies on different agents, validation loss on different tasks for agents and
more network weights visualization results.

Figure 15: Validation loss on different agents and targets. The plot shows the
validation loss on different targets combinations for agents. The target
velocity ranges from -0.08 to 0.08 and the target heights are 0.1, 0.15 and
0.20. Each subplot shows the task, running and jumping loss for one agent
given a pair of specified target velocity and height.

Figure 16: Summary of the ablation study for agents. The subplots show the normalized validation loss for task (weighted summation of different losses), running and jumping from top to bottom row. The subplots in different columns show the results of different agents. Table 3:  Ablation Summary on _Monster_. Name | Training Setting | Norm. Valid. Loss (Monster)
---|---|---
OP | AF | BS | PS | SV | TG | LD | Run. | Jump. | Task
Full | Adam | sin | 32 | ✓ | ✓ | ✓ | ✓ | 0.46$\pm$0.01 | 0.77$\pm$0.00 | 0.71$\pm$0.01
Full-OP | SGD | sin | 32 | ✓ | ✓ | ✓ | ✓ | 0.99$\pm$0.01 | 0.99$\pm$0.00 | 0.99$\pm$0.01
Full-AF | Adam | tanh | 32 | ✓ | ✓ | ✓ | ✓ | 0.65$\pm$0.01 | 0.84$\pm$0.02 | 0.79$\pm$0.02
Full-PS | Adam | sin | 32 | ✗ | ✓ | ✓ | ✓ | 0.93$\pm$0.00 | 0.99$\pm$0.00 | 0.98$\pm$0.00
Full-SV | Adam | sin | 32 | ✓ | ✗ | ✓ | ✓ | 0.51$\pm$0.01 | 0.86$\pm$0.01 | 0.78$\pm$0.01
Full-TG | Adam | sin | 32 | ✓ | ✓ | ✗ | ✓ | 0.99$\pm$0.00 | 0.97$\pm$0.01 | 0.98$\pm$0.01
Full-LD | Adam | sin | 32 | ✓ | ✓ | ✓ | ✗ | 0.99$\pm$0.00 | - | 0.90$\pm$0.06
Table 4:  Ablation Summary on _HugeStool_. Name | Training Setting | Norm. Valid. Loss (HugeStool)
---|---|---
OP | AF | BS | PS | SV | TG | LD | Run. | Jump. | Task
Full | Adam | sin | 32 | ✓ | ✓ | ✓ | ✓ | 0.40$\pm$0.01 | 0.58$\pm$0.01 | 0.53$\pm$0.01
Full-OP | SGD | sin | 32 | ✓ | ✓ | ✓ | ✓ | 0.97$\pm$0.01 | 0.96$\pm$0.01 | 0.96$\pm$0.01
Full-AF | Adam | tanh | 32 | ✓ | ✓ | ✓ | ✓ | 0.49$\pm$0.01 | 0.59$\pm$0.01 | 0.56$\pm$0.01
Full-PS | Adam | sin | 32 | ✗ | ✓ | ✓ | ✓ | 0.99$\pm$0.00 | 0.99$\pm$0.00 | 0.99$\pm$0.00
Full-SV | Adam | sin | 32 | ✓ | ✗ | ✓ | ✓ | 0.40$\pm$0.01 | 0.62$\pm$0.01 | 0.57$\pm$0.00
Full-TG | Adam | sin | 32 | ✓ | ✓ | ✗ | ✓ | 0.99$\pm$0.00 | 0.90$\pm$0.00 | 0.93$\pm$0.00
Full-LD | Adam | sin | 32 | ✓ | ✓ | ✓ | ✗ | 0.99$\pm$0.01 | - | 0.70$\pm$0.01
Table 5:  Ablation Summary on _Stool_. Name | Training Setting | Norm. Valid. Loss (Stool)
---|---|---
OP | AF | BS | PS | SV | TG | LD | Run. | Jump. | Task
Full | Adam | sin | 32 | ✓ | ✓ | ✓ | ✓ | 0.50$\pm$0.01 | 0.58$\pm$0.04 | 0.56$\pm$0.03
Full-OP | SGD | sin | 32 | ✓ | ✓ | ✓ | ✓ | 0.78$\pm$0.02 | 0.76$\pm$0.02 | 0.76$\pm$0.02
Full-AF | Adam | tanh | 32 | ✓ | ✓ | ✓ | ✓ | 0.60$\pm$0.01 | 0.63$\pm$0.03 | 0.62$\pm$0.04
Full-PS | Adam | sin | 32 | ✗ | ✓ | ✓ | ✓ | 0.98$\pm$0.01 | 0.96$\pm$0.00 | 0.96$\pm$0.02
Full-SV | Adam | sin | 32 | ✓ | ✗ | ✓ | ✓ | 0.50$\pm$0.03 | 0.64$\pm$0.05 | 0.61$\pm$0.06
Full-TG | Adam | sin | 32 | ✓ | ✓ | ✗ | ✓ | 0.98$\pm$0.00 | 0.87$\pm$0.01 | 0.89$\pm$0.01
Full-LD | Adam | sin | 32 | ✓ | ✓ | ✓ | ✗ | 0.99$\pm$0.01 | - | 0.83$\pm$0.05
Ouput Hidden | sin | tanh | relu | gelu | sigmoid
---|---|---|---|---|---
sin | 0.745 | 0.771 | 0.765 | 0.803 | 0.810
tanh | 0.753 | 0.805 | 0.783 | 0.799 | 0.849
relu | 0.780 | 0.770 | 0.806 | 0.775 | 0.889
gelu | 0.767 | 0.780 | 0.773 | 0.795 | 0.829
sigmoid | 0.754 | 0.767 | 0.800 | 0.799 | 0.860
Table 6:  Full ablation studies for activation functions. The values in the
table represent the normalized validation loss for _Task_. For each setting,
the experiments are repeated for multiple times. It can be observed that the
model using $sin$ for both hidden and output layer achieves the best results.


#  Numerical bifurcation analysis of improved denatured Morris-Lecar neuron
model

H.O. Fatoyinboa*, S.S. Munia, I. Ghosha, I.O. Sarumib, A. Abidemic

aSchool of Mathematical and Computational Sciences
Massey University, Palmerston North, New Zealand
bDepartment of Mathematics and Statistics
King Fahd University of Petroleum and Minerals, Dhahran, Saudi Arabia
cDepartment of Mathematical Sciences
Federal University of Technology, Akure, Nigeria

###### Abstract

It is well-known that the electrical activities of neurons are induced by a
wide variety of external factors. This work considers the effect of
electromagnetic induction on improved denatured Morris-Lecar neuron model. The
dependence of dynamical behaviour of the original denatured Morris-Lecar model
on parameters is addressed through numerical bifurcation analysis. This allows
us to explore the changes in dynamics of the model qualitatively as parameters
are varied. Then we investigate the effects of external periodic current and
electromagnetic flux on the dynamical properties of the improved denatured
Morris-Lecar neuron model. Different types of dynamical behaviour, ranging
from regular periodic spiking to complex bursting, are found when the multiple
parameters are varied simultaneously. The improved model could be applied to
research where simple models are required for physiological and
pathophysiological responses in neurons.

††E-mail addresses: H. O. Fatoyinbo (h.fatoyinbo@massey.ac.nz), S.S. Muni
(s.muni@massey.ac.nz), I. Ghosh(i.ghosh@massey.ac.nz), I. O.
Sarumi(isarumi@kfupm.edu.sa), A. Abidemi<EMAIL_ADDRESS>
##  1 Introduction

Nonlinear dynamical behaviours of neurons have been widely studied to
understand the mechanisms associated with neuronal diseases such as epilepsy
[1] and Parkinson’s disease [2]. They are frequently modelled by a nonlinear
system of differential equations. In 1952, Hodgkin and Huxley (H-H)
constructed a neuronal model from experimental studies to describe the
conduction and propagation of electrical impulses in the squid giant axon [3].
Following the H-H model, FitzHugh and Nagumo [4], [5] constructed a two-
variable model which is a simplified version of the H-H model to describe the
different phenomena in neuron dynamics. Other well-known models are the Morris
and Lecar (ML) model [6], the Hindmarsh–Rose model [7], and the Izhikevich
neuron model [8].

Bifurcation analysis has been employed in many previous studies to uncover
transitions of firing patterns in neurons and other biological cells [11],
[28], [29]. It is very helpful to identify important parameters that can cause
significant change in the dynamical behaviour qualitatively. Firing activities
of neurons are classified into different classes of excitability depending on
the transitions from a rest state to a firing (oscillatory) state. The authors
in [13] classified neurons into two classes of excitability based on the
onsets of firing. In class I excitability, the periodic oscillations occur
through a saddle-node on an invariant circle (SNIC) bifurcation. In contrast,
for class II excitability the periodic oscillations emerge through a Hopf
bifurcation. The mode of transitions between the two classes of excitability
is characterised by bifurcations. On this basis, the codimension-two
bifurcations associated with a change between classes of excitability are
identified in previous studies of neurons and other excitable cells [30],
[31], [16], [32].

Multiple researchers have studied various factors that influence the firing
activities in neurons and other excitable cells [9],[10],[11],[12]. The effect
of direct current stimulus [14], temperature [15], pressure [16], time-delay
[17] and environmental noise [18] have been investigated. Magnetic induction
is another factor that affects the neuron membrane potential significantly,
which has been observed in physiological experiments [19],[20]. According to
Maxwell’s electromagnetic induction theorem, variation of neuronal membrane
potential results in the generation of a magnetic field in the neuron and the
induction current can influence the firing activities of the neuron. Li et al.
[21] constructed a mathematical model to study the variation of firing
patterns of neurons under the influence of electromagnetic radiation. Lv and
Ma [22] used the magnetic flux to describe the effects of electromagnetic
radiation on the electrical activities of neurons by coupling the magnetic
flux with the membrane potential. Many recent papers have studied the effects
of electromagnetic induction on neurons and collective behaviours in a network
of neurons [24], [25], [26], [27]. Diverse firing modes and complicated
dynamical behaviour including bursting are observed in neuron models under
electromagnetic radiation. Also, spatiotemporal dynamics including chaotic
bursting, pattern formation and chimera state are found in a network of
neurons [23].

In many studies, high-dimensional conductance-based models have been used to
explore electrophysiological activities in neurons. However, due to the
nonlinear behaviour of the models, they are computationally expensive; hence
low-dimensional models have been instrumental to study electrical activities
in neurons because the reduced models capture the excitable properties in
higher-dimensional models. A two-variable ML model is a reduced model that
replicates the dynamical behaviours observed in the 8-dimensional H-H model.
In [33], a simplified version of the ML model, known as the denatured ML
neuron model, is introduced. To the best of our knowledge, this is the first
work to give a more comprehensive analysis of the denatured ML model

In this paper, we will use the denatured ML model to investigate the dynamics
of neurons under the influence of electromagnetic induction and external
periodic current. The paper is structured as follows: In Section 2, we
describe the model and the parameters. In Section 3, we analyse the original
denatured ML model via phase plane and bifurcation analyses to explore the
firing patterns of neurons. We also explore the effects of electromagnetic
flux and external periodic current on firing patterns of the improved ML model
in Section 4. The summary and conclusion of results are given in Section 5.

##  2 Model Description

The denatured Morris-Lecar neuron model is a simplified version of the two
variable Morris-Lecar neuron model [6]. It consists of two-coupled ODEs
described in [33] as

$\displaystyle\dot{x}$ $\displaystyle\\!\begin{aligned}
&=x^{2}(1-x)-y+I,\end{aligned}$ (2.1a) $\displaystyle\dot{y}$
$\displaystyle\\!\begin{aligned} &=Ae^{\alpha x}-\gamma y,\end{aligned}$
(2.1b)

where $x$ is the membrane action potential and $y$ is the recovery variable.
The system parameter $I$ is the external current, $A$, $\alpha$ and $\gamma$
are positive constants. In this paper we improve model (2.1) by adding
electromagnetic induction and external periodic current. The dynamical
equations of the improved version of (2.1) are described as follows

$\displaystyle\dot{x}$ $\displaystyle\\!\begin{aligned}
&=x^{2}(1-x)-y+I+k\rho(\phi)x,\end{aligned}$ (2.2a) $\displaystyle\dot{y}$
$\displaystyle\\!\begin{aligned} &=Ae^{\alpha x}-\gamma y,\end{aligned}$
(2.2b) $\displaystyle\dot{\phi}$ $\displaystyle\\!\begin{aligned}
&=k_{1}x-k_{2}\phi+\phi_{\rm ext},\end{aligned}$ (2.2c)

where $k$ is the feedback gain, the function
$\rho(\phi)=\alpha_{1}+3\beta\phi^{2}$ is the electromagnetic effect on the
neuron and $\phi$ is the magnetic flux across the cell membrane. $\alpha_{1}$
and $\beta$ are the memory conductances, and $\phi_{\rm ext}$ is the external
magnetic flux. The external periodic current $I=I_{0}\sin(\omega t)$, where
$I_{0}$ is the amplitude of the current and $\omega$ is the angular frequency.
Unless otherwise stated, the parameter values except for $I_{0}$ and $\gamma$
are fixed as: $A=0.0041$, $\alpha=5.276$, $k=0.003$, $k_{1}=0.19$,
$k_{2}=0.5$, $\alpha_{1}=0.1$, $\beta=0.02$, $\omega=0.01$, and $\phi_{\rm
ext}=0.2$.

##  3 Analysis of the original denatured Morris-Lecar model

We begin our analysis by investigating the dynamical behaviour of the original
denatured ML system (2.1) via the interaction of the variables in the phase
plane. The $x$\- and $y$-nullclines are set of points in the phase plane such
that $\dot{x}=0$ and $\dot{y}=0$, that is $x$ and $y$ are constant, and the
intersection of the nullclines defines the equilibrium point. The
$y$-nullcline is shown in black and the $x$-nullclines are shown in orange for
various values of $I$ in Fig. 1. The blue [red] filled circle corresponds to
stable [unstable] equilibrium. For extremely low values of $I$ the two
nullclines intersect at one point, corresponding to a unique stable
equilibrium. For example, the dashed orange curve in Fig. 1 is when
$I=-0.025$. For intermediate values of $I$ the nullclines intersect at three
points, corresponding to three equilibria of the system. The case of $I=0.013$
is the solid orange curve in Fig. 1, there is one stable and two unstable
equilibria. For large values of $I$ the two nullclines intersect at one point,
this corresponds to a unique equilibrium point. The point is either stable or
unstable depending on the value of $I$. For $I=0.07$, the equilibrium point is
stable, see the dash-dot orange curve in Fig. 1. From the analysis above, it
is observed that the number and stability of equilibria of the system change
as the external current $I$ is varied, thus to further our understanding of
how the model dynamics changes as parameters are varied we will consider
numerical bifurcation analysis of system (2.1) in the next section.

Figure 1: The $(x,y)$ phase plane of (2.1). The $y$-nullcline is plotted black
and the $x$-nullclines are plotted orange for three different values of $I$.
The dashed, solid and dashed-dot curves are $x$-nullclines for $I=-0.025$,
$0.013$, and $0.07$, respectively. The blue [red] filled circle correspond to
stable and unstable equilibria.

###  3.1 One-parameter bifurcation analysis

In this section, we consider the dependence of the dynamical behaviour of
(2.1) on parameters $I$ and $\gamma$ through numerical bifurcation analysis.
The bifurcation diagrams were produced in MatCont, a numerical continuation
software [34]. Fig. 2(a) is a one-parameter bifurcation diagram of the model
(2.1) for the bifurcation parameter $I$. The blue [red] curves correspond to
stable [unstable] solutions, respectively. The points labelled ${\rm LP}_{1}$
and ${\rm LP}_{2}$ are saddle-node bifurcation points. The points labelled
${\rm HB}_{1}$ and ${\rm HB}_{2}$ are Hopf bifurcation points. For relatively
low values of $I$ there exists a unique stable solution. With increasing $I$,
the stable solution loses stability in a saddle-node bifurcation ${\rm
LP}_{1}$ at $I=0.0153$. Upon further increase in the value of $I$, the
unstable solution regains stability in a saddle-node bifurcation ${\rm
LP}_{2}$ at $I=0.0109$. Between the two saddle-node bifurcations, ${\rm
LP}_{1}$ and ${\rm LP}_{2}$, there exist three solutions: two stable and one
unstable. Thus, in this parameter regime there is bistability (two stable
solutions coexist) in the system. The stable solution produced at ${\rm
LP}_{2}$ loses stability in a Hopf bifurcation ${\rm HB}_{1}$ at $I=0.0117$.
Between saddle-node ${\rm LP}_{1}$ and Hopf bifurcation ${\rm HB}_{1}$ there
are three solutions: one stable (lower branch) and two unstable (middle and
upper branch). An example of this scenario is shown in Fig. 1 when $I=0.013$.
Since the first Lyapunov coefficient of ${\rm HB}_{1}$ is negative
($a=-5.2312$) then ${\rm HB}_{1}$ is a supercritical bifurcation. As known
from the theory of supercritical Hopf bifurcation [35], a stable periodic
solution emanates from ${\rm HB}_{1}$. As the value of $I$ increases the
stable periodic solution terminates in a Hopf bifurcation ${\rm HB}_{2}$ at
$I=0.0541$. The first Lyapunov coefficient of ${\rm HB}_{2}$ is negative
($a=-0.2643$), thus it is a supercritical Hopf bifurcation. Also, the unstable
solution from ${\rm HB}_{1}$ regains stability at ${\rm HB}_{2}$. For
relatively high values of $I$ there is a unique stable solution.

(a)

(b)

Figure 2: Bifurcation diagrams of (2.1) with (a) $I$; (b) $\gamma$ as
bifurcation parameters. Blue [red] curves correspond to stable [unstable]
solution. HB: Hopf bifurcation; LP: saddle-node bifurcation (of an
equilibrium).

Next $\gamma$ is varied. A bifurcation diagram is shown in Fig. 2(b). We
observed that for the parameter values of $\gamma$ considered, the system has
a unique solution. Apart from the parameter regime between the two Hopf
bifurcations ${\rm HB}_{1}$ and ${\rm HB}_{2}$ where the solution is unstable,
it is stable everywhere else (to the left ${\rm HB}_{1}$ and the right of
${\rm HB}_{2}$). The first Lyapunov coefficients of ${\rm HB}_{1}$ and ${\rm
HB}_{2}$ are negative hence they are supercritical Hopf bifurcations. The
stable periodic solution that emanates from ${\rm HB}_{1}$ terminates at ${\rm
HB}_{2}$. From the above analysis, system (2.1) models neuron activity with
class II excitability.

###  3.2 Two-parameter bifurcation analysis

Here a two-parameter bifurcation analysis of model (2.1) by varying $I$ and
$\gamma$ simultaneously is considered. Fig. 3 is a codimension-2 bifurcation
diagram of (2.1) in $(I,\gamma)$-plane. The black curves correspond to the
loci of the Hopf and saddle-node bifurcations in Fig. 2(a). The saddle-node
loci collide and annihilate in a codimension-2 bifurcation known as cusp
bifurcation CP. A generalised Hopf bifurcation GH appears along the Hopf
locus, this is a codimesion-2 point where the Hopf bifurcation changes from
subcritical to supercritical. Fig. 3 is divided into four qualitatively
different dynamical behaviour regions, with enlargement in Fig. 3b. Each
region is assigned a number and a colour. The remainder of this section
describes the dynamical behaviour in each region.

Figure 3: A two-parameter bifurcation diagram of model (2.1) in the
$(I,\gamma)$-plane; (b) is the enlargement of Fig. 3a. The black curves are
the loci of the saddle-node and Hopf bifurcations in Fig. 2. The codimension-2
bifurcation points are CP: cusp bifurcation; GH: generalised Hopf bifucation.

In region \scriptsize{1}⃝, the system has one stable equilibrium and no
periodic solutions (rest state). In region \scriptsize{2}⃝, the system has
three equilibria: two stable and one unstable equilibria. This region is
bounded by the loci of saddle-node bifurcations, thus there is no periodic
solution. Apart from the loci of the saddle-node bifurcations previously
observed there is now a Hopf bifurcation locus between which there exists a
stable periodic solution in region \scriptsize{3}⃝. Finally, in region
\scriptsize{4}⃝ there is one unstable equilibrium and one stable periodic
solution.

##  4 Analysis of improved denatured Morris-Lecar model

In this section, the dynamics of the improved denatured ML model (2.2) under
the influence of electromagnetic induction and external periodic current is
investigated . The numerical simulations were carried out in Python using
odeint solver, which uses a variable step-size Runge-Kutta method. The time
evolution of the membrane potential $x$ and the projection of the phase space
in $(x,y)$-plane are shown in Fig. 4. Various firing activities including
regular spiking and bursting are found as the amplitude of the external
current $I_{0}$ is varied and set $\gamma=0.315$. For $I_{0}=0.00072$, the
model is in a rest state (see Fig. 4(a)). With the increase of $I_{0}$,
regular periodic spiking activity is observed, an example is shown in Fig.
4(c) for $I_{0}=0.0155$. However, increasing $I_{0}$ further the model shows
regular periodic bursting. The model has period-2, period-4 and period-13 when
$I_{0}=0.016,0.02,$ and $0.04$, respectively (see Figs. 4e, g and i). As seen
from the above analysis, complex dynamical behaviours occur in the improved
denatured ML model under the influence of the amplitude of the external
current $I_{0}$.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

Figure 4: Time series plots of $x$ for model (2.2) for some values of $I_{0}$
with $\gamma=0.315$. (a) $I_{0}=0.00072$; (c) $I_{0}=0.0155$; (e)
$I_{0}=0.016$; (g) $I_{0}=0.02$; (i) $I_{0}=0.04$. Their corresponding phase
plane in $(x,y)$-plane are shown in (b), (d), (f), (h), and (j), respectively.

Next, multiple parameters are varied simultaneously to explore the effects of
the electromagnetic radiation on the dynamical behaviors of the improved
denatured ML model. As seen in Fig. 5, model (2.2) exhibits mixed mode
oscillations and bursting activities. In Fig. 5a-b, we observe mixed mode
oscillations where several regular single spikes and subthreshold oscillations
coexist. Both of these cases have $A=0.005$ with variations in
$\gamma,\omega,k$ and $I_{0}$. For
$A=0.002,\gamma=0.1576,\omega=0.02,k=0.018$, and $I_{0}=0.17$, there exists
period-7 bursting with subthreshold oscillations (see Fig. 5c). In Fig. 5d, we
observe that when $k=0.919$, $\omega=0.001$ with $A=0.002$, $\gamma=0.1576$ ,
and $I_{0}=0.17$ the dynamics subsides to periodic spiking. Finally in Fig.
5e-f regular periodic bursting patterns appear.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 5: Time series plots of $x$ for model (2.2) for (a)
$A=0.005,\gamma=0.23,\omega=0.03$, and $I_{0}=0.04$; (b)
$A=0.005,\gamma=0.035,\omega=0.0121,k=0.0231$, and $I_{0}=0.187$; (c)
$A=0.002,\gamma=0.1576,\omega=0.02,k=0.018$, and $I_{0}=0.17$; (d)
$A=0.002,\gamma=0.1576,\omega=0.001,k=0.919$, and $I_{0}=0.17$; (e)
$A=0.0187,\gamma=0.231,\omega=0.015,k=0.05$, and $I_{0}=0.434$; (f)
$A=0.0018,\gamma=0.231,\omega=0.015,k=0.08$, and $I_{0}=0.201$.

##  5 Conclusion

In this study, phase plane analysis and bifurcation analysis are employed to
study the firing activities in a simple electrophysiological model of neurons.
By using codimension-1 and -2 bifurcation analysis, we showed how variation of
model parameters affects the dynamics of the membrane potential such as
transitions from a rest state to periodic oscillations and vice versa. We
found that model exhibits class II excitability. The class of excitability
affects how neurons decode and process information, thus it is important to
understand the mechanisms underlying the transitions between classes of
excitability. As in [33], the denatured ML model exhibits class I
excitability. In future work we will investigate the mechanisms that govern
switches between classes of excitability in the improved model.

In addition, we investigated the dynamics of the improved model behaviour
under the influence of external periodic current and electromagnetic
induction. Evidence from numerical simulations shows that the electrical
activity of the improved neuron model depends greatly on system parameters. We
observed that for different values of the external periodic current the model
produces simple and complex oscillations including multiple mode bursting
activities. Also, we explored the effects of the angular frequency $\omega$
and the feedback gain $k$ on the behaviour of the improved model.

The model displays qualitative behaviours observed in experiments and
computational analysis of higher dimension neuron models. This work provide
for some applications where simple models are required for physiological and
pathophysiological responses in neurons and other excitable cells.

In the present study, we have considered analysis of a single neuron, however,
neurons are interconnected through chemical and electrical synapes. Hence, the
analysis of ring network of neurons under electromagnetic induction will be a
subject of future work.

## References

  * [1] G. F. Ayala, M. R. Dichter, J. Gumnit, H. Matsumoto, and W. A. Spencer, “Genesis of epileptic interictal spikes. New knowledge of cortical feedback systems suggests a neurophysiological explanation of brief paroxysms,” _Brain Res_ , vol. 52, p. 1–17, 1973.
  * [2] R. Levy, W. Hutchison, A. Lozano, and J. Dostrovsky, “High-frequency synchronization of neuronal activity in the subthalamic nucleus of Parkinsonian Patients with Limb Tremor,” _J. Neurosci_ , vol. 20, p. 7766–7775, 2000.
  * [3] A. L. Hodgkin, “The local electric changes associated with repetitive action in a non-medullated axon,” _J Physiol_ , vol. 107, no. 2, pp. 165–181, 3 1948.
  * [4] R. FitzHugh, “Impulses and physiological states in theoretical model of nerve membrane,” _Biophysical J_ , vol. 1, pp. 445–466, 1961.
  * [5] J. Nagumo, S. Arimoto, and S. Yoshizawa, “An active pulse transmission line simulating nerve axon,” _Proceedings of the IRE_ , vol. 50, no. 10, pp. 2061–2070, 1962.
  * [6] C. Morris and H. Lecar, “Voltage oscillations in the barnacle giant muscle fiber,” _Biophysical J_ , vol. 35, pp. 193–213, 1981.
  * [7] J. L. Hindmarsh and R. M. Rose, “A model of neuronal bursting using three coupled first order differential equations,” _Proc. R . Soc. Lond. B_ , vol. 221, pp. 87–102, 1984.
  * [8] E. M. Izhikevich, _Dynamical systems in neuroscience : the geometry of excitability and bursting_. Cambridge: MIT Press, 2007.
  * [9] Z. Zhao and H. Gu, “Transitions between classes of neuronal excitability and bifurcations induced by autapse,” _Scientific Reports_ , vol. 7, no. 1, p. 6760, 7 2017.
  * [10] C. Liu, X. Liu, and S. Liu, “Bifurcation analysis of a Morris–Lecar neuron model,” _Biol Cybern_ , vol. 108, no. 1, pp. 75–84, 2 2014.
  * [11] W. Govaerts and B. Sautois, “The onset and extinction of neural spiking: A numerical bifurcation approach,” _J Comput Neurosci_ , vol. 18, no. 3, pp. 265–274, 6 2005.
  * [12] J. M. González-Miranda, “Pacemaker dynamics in the full Morris-Lecar model,” _Commun Nonlinear Sci Numer Simul_ , vol. 19, pp. 3229–3241, 2014\.
  * [13] J. Rinzel and G. B. Ermentrout, “Analysis of neural excitability and oscillations,” in _C. Koch, I. Segev 2nd (Eds) , Methods in neuronal modeling: From ions to network_ , C. Koch and I. Segev, Eds. London: MIT Press, 1999, pp. 251–291.
  * [14] H. Wang, L. Wang, L. Yu, and Y. Chen, “Response of Morris-Lecar neurons to various stimuli,” _Phys. Rev. E_ , vol. 83, p. 021915, 2011.
  * [15] M. Xing, X. Song, Z. Yang, and Y. Chen, “ Bifurcations and excitability in the temperature-sensitive Morris–Lecar neuron,” _Nonlinear Dyn_ , vol. 100, p. 2687–2698, 2020.
  * [16] H. O. Fatoyinbo, R. G. Brown, D. J. W. Simpson, and B. van Brunt, “Numerical bifurcation analysis of pacemaker dynamics in a model of smooth muscle cells,” _Bull Math Bio_ , vol. 82, no. 95, 2020.
  * [17] B. Jia, Y. Wu, D. He, B. Guo, and L. Xue, “Dynamics of transitions from anti-phase to multiple in-phase synchronizations in inhibitory coupled bursting neurons,” _Nonlinear Dyn_ , vol. 93, no. 9, p. 1599–1618, 2018\.
  * [18] Y. Gong, Y. Xie, X. Lin, and Y. Hao, “Non-gaussian noise‐ optimized intracellular cytosolic calcium oscillations,” _Biosystems_ , vol. 103, no. 1, p. 13–17, 2011.
  * [19] J. F. Barry, M. J. Turner, J. M. Schloss, D. R. Glenn, Y. Song, M. D. Lukin, H. Park, and R. L. Walsworth, “Optical magnetic detection of single-neuron action potentials using quantum defects in diamond,” _Proc. Natl. Acad. Sci. U.S.A._ , vol. 113, p. 14133–14138, 2016.
  * [20] A. S. Etémé and A. M. Tabi, “Firing and synchronization modes in neural network under magnetic stimulation.” _Commun Nonlinear Sci_ , vol. 72, p. 420–440, 2019.
  * [21] J. Li, Y. Wu, M. Du, and W. Liu, “Dynamic behavior in firing rhythm transitions of neurons under electromagnetic radiation,” _Acta Phys. Sin_ , vol. 64, p. 030503, 2015.
  * [22] M. Lv and J. Ma, “Multiple modes of electrical activities in a new neuron model under electromagnetic radiation,” _Neurocomputing_ , vol. 205, p. 375–381, 2016.
  * [23] A. Karthikeyan, M. I., K. Rajagopal, and P. Duraisamy, “Effect of temperature sensitive ion channels on the single and multilayer network behavior of an excitable media with electromagnetic induction,” _Chaos, Solitons & Fractals_, vol. 150, p. 111144, 2021.
  * [24] X. Yang, G. Zhang, X. Li, and D. Wang, “The synchronization behaviors of coupled fractional-order neuronal networks under electromagnetic radiation,” _Symmetry_ , vol. 13, p. 2204, 2021.
  * [25] S. S. Muni, H. O. Fatoyinbo, and I. Ghosh, “Dynamical effects of electromagnetic flux on chialvo neuron map: nodal and network behaviors,” 2022\.
  * [26] A. Mondal, R. K. Upadhyay, J. Ma, B. K. Yadav, and S. K. Sharma, “Bifurcation analysis and diverse firing activities of a modified excitable neuron model,” _Cogn Neurodyn_ , vol. 13, no. 4, pp. 393–407, 8 2019.
  * [27] K. Rajagopal, M. I., B. Ramakrishnan, A. Karthikeyan, and P. Duraisamy, “Modified morris-lecar neuron model: effects of very low frequency electric fields and of magnetic fields on the local and network dynamics of an excitable media,” _Nonlinear Dyn_ , vol. 104, p. 4427–4443, 2021.
  * [28] L. Duan and Q. Lu, “Codimension-two bifurcation analysis on firing activities in Chay neuron model,” _Chaos, Solitons & Fractals_, vol. 30, no. 5, pp. 1172–1179, 12 2006.
  * [29] H. O. Fatoyinbo, S. S. Muni, and A. Abidemi, “Influence of sodium inward current on the dynamical behaviour of modified morris-lecar model,” _Eur. Phys. J. B_ , vol. 95, no. 1, pp. 1–15, 2022.
  * [30] K. Tsumoto, H. Kitajima, T. Yoshinaga, K. Aihara, and H. Kawakami, “Bifurcations in Morris-Lecar neuron model,” _Neurocomputing_ , vol. 69, no. 4-6, pp. 293–316, 1 2006.
  * [31] L. Duan, Q. Lu, and Q. Wang, “Two-parameter bifurcation analysis of firing activities in the Chay neuronal model,” _Neurocomputing_ , vol. 72, no. 1-3, pp. 341–351, 12 2008.
  * [32] F. Wu, H. Gu, and Y. Jia, “Bifurcations underlying different excitability transitions modulated by excitatory and inhibitory memristor and chemical autapses,” _Chaos, Solitons & Fractals_, vol. 153, p. 111611, 2021.
  * [33] D. G. Schaeffer and J. W. Cain, _Ordinary Differential Equations: Basics and Beyond_. New York: Springer, 2016\.
  * [34] W. Govaerts, Y. A. Kuznetsov, V. De Witte, A. Dhooge, H. G. E. Meijer, W. Mestrom, A. M. Riet, and B. Sautois, _MATCONT and CL MATCONT: Continuation toolboxes in matlab_. Gent University and UtrechtUniversity, 2011.
  * [35] Kuznetsov, Y. A., _Elements of applied bifurcation theory_ , 3rd ed. New York: Springer-Verlag, 1995\.