# Solitary states in complex networks: impact of topology

Leonhard Schülen  Institut für Theoretische Physik, Technische Universität
Berlin, Hardenbergstr. 36, 10623 Berlin, Germany  Maria Mikhailenko  Charité –
Universitätsmedizin Berlin, Einstein Center for Neurosciences Berlin,
Charitéplatz 1, 10117, Berlin, Germany  Everton S. Medeiros  Institute for
Chemistry and Biology of the Marine Environment, Carl-von-Ossietzky-University
Oldenburg, Carl-von-Ossietzky-Straße 9 - 11, 26111 Oldenburg, Germany  Anna
Zakharova  Institut für Theoretische Physik, Technische Universität Berlin,
Hardenbergstr. 36, 10623 Berlin, Germany  Bernstein Center for Computational
Neuroscience, Humboldt-Universität zu Berlin, Philippstraße 13, 10115 Berlin,
Germany

###### Abstract

The dynamical behavior of networked systems is expected to reflect the
features of their coupling structure. Yet, symmetry-broken solutions often
occur in symmetrically coupled networks. An example is provided by the so-
called solitary states where the dynamics of one network node is different
from the entire symmetric network. Here, we investigate the structural
constraints of networks for the appearance of solitary states in their
dynamics. By performing a large number of numerical simulations, we find that
such states occur with high probability in asymmetric networks, such as the
ones exhibiting the scale-free property. Next, we analyze the structural
features of the networks demonstrating solitary states to reveal that the
minimum number of connections in the adjacent nodes of a solitary one is
crucial for the appearance of the solitary states. Finally, we perform
bifurcation analysis of a dimension-reduced system, confirming the influence
of the connectivity of the neighboring nodes on the solitary ones.

##  1 Introduction

Understanding the interplay between network topology and emerging dynamics is
one of the central issues in the field of nonlinear dynamics and the theory of
complex networks [1, 2, 3]. Often, perfectly symmetric networked systems
exhibit symmetry-broken solutions [4, 5]. For networks of oscillators, a
particularly interesting example is the extreme case of cluster
synchronization in which only one oscillator does not synchronize with the
rest of the network. This highly unbalanced network configuration is called
“solitary” state [6, 7]. The occurrence of these states has been observed in a
variety of symmetrically coupled systems such as globally coupled nonlinear
oscillators [6], in non-locally (and globally) coupled networks of phase
oscillators [8, 9], coupled chaotic maps [10, 11, 12], coupled excitable
systems [13], multilayer networks [14, 15, 16, 17], and time-delayed systems
[18]. Such ubiquity naturally drew attention to solitary states and raised
questions about the mechanism behind their appearance [19, 15]. Recently,
Schülen et al. have addressed this issue for coupled neural oscillators and
found that solitary states are created sub-critically in a fold-bifurcation
[20]. Additionally, they have demonstrated that the dynamics of the solitary
oscillator can also be chaotic while the synchronized cluster stays periodic
[20].

Despite the advanced understanding of ubiquity and onset of the solitary
states, some essential questions remain unanswered. For instance, what are the
topological constraints for the emergence of these asymmetric solutions?
Naturally, symmetry-broken solutions are unexpected in the scope of symmetric
systems with uniform distribution of node degree. However, complex networks
with asymmetries in their structure of connections, such as the ones with
scale-free property presenting a power-law distribution of node degree, do not
offer any immediate advantage for the occurrence of solitary states. In fact,
it has been recently found that asymmetries can actually favor synchronization
[21, 22] or regular behavior [23, 24]. Therefore, the onset of solitary states
in asymmetric networks is as equally important as the symmetric case,
especially for a better understanding of how these states occur in general
setups.

In this work, we address the onset, and the topological dependencies of
solitary states in complex networks of FitzHugh-Nagumo oscillators. We begin
our investigation by studying the appearance of solitary states in three
different network topologies: a symmetric nonlocally coupled ring, a random,
and a scale-free graph. To establish the range of parameters in which the
solitary states would appear, we use the symmetric nonlocally coupled ring as
a baseline for the parameter search. We find that the random and the scale-
free networks exhibit solitary states for similar parameter values as the
symmetric one. These results are verified over an ensemble of realizations of
such networks with different initial conditions. Further, we proceed to the
topological analysis of the networks, aiming to identify a common feature of
solitary nodes. We analyze their node degrees, neighbor node degrees, and
eigenvector centrality. Among these measures, we found that a high neighbor
node degree correlates with the onset of solitary states. Our next step is to
verify this observation more rigorously by performing bifurcation analysis.
For that, we consider geometrical arguments to obtain a dimension-reduced
system, where the bifurcation analysis demonstrates the onset of the solitary
states in a fold bifurcation and their dependencies on the minimum value of
the neighbors node degree.

##  2 Solitary states in complex networks

We consider different network topologies, each network consisting of identical
FitzHugh-Nagumo (FHN) elements in the oscillatory regime. Such networks are
numerically simulated using the following dimensionless equations:

$\displaystyle\varepsilon\frac{du_{i}}{dt}$
$\displaystyle=f(u_{i},v_{i})+\sum\limits_{j=1}^{N}\frac{A_{ij}}{d_{i}}[\sigma_{u}(u_{j}-u_{i})+\sigma_{v}(v_{j}-v_{i})],$
$\displaystyle\frac{dv_{i}}{dt}$ $\displaystyle=g(u_{i},v_{i}),$ (1)

where the pair ($u_{i}$, $v_{i}$) accounts for the activator and inhibitor
variables, respectively, of each FHN oscillator $i$ with $i=1,\dots,N$. The
parameter $N$ prescribes the network size. The functions:

$f(u_{i},v_{i})=u_{i}-\frac{u_{i}^{3}}{3}-v_{i},\quad
g(u_{i},v_{i})=g(u_{i})=u_{i}+a$ (2)

specify the local dynamics of the variables $u_{i}$ and $v_{i}$. The parameter
$\varepsilon$ defines the time scale separation between the fast ($u_{i}$) and
slow ($v_{i}$) variables. Throughout this study, we fix $\varepsilon=0.1$. For
isolated oscillators, the threshold parameter $a$ separates the excitatory
($|a|>1$) from the oscillatory ($|a|<1$) regime through a supercritical Hopf
bifurcation. We set each oscillator in the oscillatory regime by keeping
$a=0.5$. The different network topologies are characterized by their adjacency
matrix $A_{ij}$. The FHN oscillators are diffusively coupled via both
variables $u_{i}$ and $v_{i}$, similarly to [19, 25, 26, 27, 28]. The strength
of the diffusive coupling is given by $\sigma_{u}$ and $\sigma_{v}$ for the
activator and the inhibitor variable, respectively. The coupling function of a
given FHN oscillator $i$ is normalized by the number of oscillators coupled to
it, i.e., its node degree $d_{i}=\sum_{j=1}^{N}A_{ij}$.

In our study, we consider three different network topologies, namely, a
symmetric nonlocally coupled ring, a random, and a scale-free network. The
topology of the symmetrically coupled ring is determined by a coupling radius
$R$ common to all nodes in the network. The coupling radius $R$ defines the
number of neighbors in each direction on a ring. This coupling scheme yields a
uniform node degree $d_{i}=2R$ $\forall$ $i\in[1,N]$, see the corresponding
graph in Fig. 1(a). This symmetric case is used as a reference case for the
comparison with the asymmetric structures. To obtain the random topology, we
consider the well-known Erdős-Rényi (ER) algorithm [29]. Starting from $N=100$
unconnected nodes, we establish a connection between any given pair of nodes
with a probability $p=0.18$. This procedure results in a normal distribution
of node degree with average value of $\langle d_{i}\rangle=18$, see the
corresponding graph in Fig. 1(b). Next, in order to obtain a scale-free
network, we use the Barabási-Albert (BA) algorithm which considers
preferential attachment in the network growth process [30, 31]. Starting from
a network with a star-like topology, each step of the BA algorithm consists of
adding $m$ new links into such a network. The newly added links are
distributed across the network as follows: a given node $i$ receives a new
link in accordance with its current degree $d_{i}$ with the following
probability $p_{i}=\frac{d_{i}}{\sum_{i}d_{i}}$. This procedure yields a
network with a power-law distribution of node degree. In Fig. 1(c), we show
the resulting graph of the BA algorithm with $m=10$.

Figure 1: Example of a solitary state in: (a) a nonlocally coupled ring graph
(for better visibility we use $N=30$ and $R=5$, yielding $d_{i}=d=10$); (b) an
Erdős-Rényi graph with $N=100$, $p=0.18$, yielding an average node degree of
$\langle d_{i}\rangle\approx 18$; (c) Barabási-Albert graph with $N=100$ and
$m=5$ links added at each algorithm step, yielding $\langle
d_{i}\rangle\approx 9.5$. The size of the nodes is proportional to the
respective node degree.

We now investigate the onset of solitary states in the considered
asymmetrically coupled networks. For that, we attribute random initial
conditions (ICs) to the system in Eq. (1) with the adjacency matrix $A_{ij}$
first prescribing the scale-free topology shown in Fig. 1(c). With this, in
Fig. 2(a), we show the time evolution of the network, illustrating the
presence of one solitary node split off from the main synchronized cluster.
This behavior can be further visualized in a snapshot of the variable $u_{i}$
at $t=2000$ (arb. units) on display in Fig. 2(b). In general, the trajectories
of solitary nodes approach an attractor coexisting with the one hosting the
synchronized cluster. The attractor of a solitary node can be periodic or
chaotic [19, 20]. For the case depicted in Figs. 2(a) and 2(b), the dynamics
of the solitary node follows a period-$3$ limit cycle as shown in the state-
space projection ($u_{i}$,$v_{i}$) depicted in Fig. 2(c). This figure
demonstrates the existence of solitary states in asymmetric networks.

Figure 2: Solitary state with a single solitary node in a scale-free network
with $N=100$, $m=5$ (as shown in Fig. 1c), $\sigma_{u}=0.12$, and
$\sigma_{v}=0.15$: (a) Space-time plot for the activator variable $u_{i}$; (b)
Snapshot of the activator variable $u_{i}$ taken at $t=2000$ (arb. units).
Blue circles correspond to the synchronized cluster, while the red one shows
the solitary node; (c) State-space projections for the solitary node (red) and
the synchronized cluster (blue).

In order to establish the coupling intensities for which solitary states arise
in complex networks, we obtain the map of regimes in the
($\sigma_{u},\sigma_{v}$)-parameter plane for all three topologies shown in
Fig. 1. For that, we consider an ensemble of $100$ network realizations with
different ICs. For each simulation, we analyze a time interval of $\Delta
t=2000$ (arb. units) for establishing the occurrence of solitary states. This
approach provides statistical results in which the fraction of simulations
resulting in solitary states is characterized by the probability of their
occurrence. With this, in the map of regimes shown in Figs. 3(a)-(c) for the
respective topologies, the color-code stands for the probability of obtaining
a solitary state for a given parameter pair ($\sigma_{u},\sigma_{v}$). We
observe that the parameter regions with a high probability of finding solitary
states occur for very similar values of the coupling intensities $\sigma_{u}$
and $\sigma_{v}$ for all three topologies. In addition, the shape of such
parameter regions is similar across the topologies. Therefore, the analysis of
the ($\sigma_{u},\sigma_{v}$)-parameter plane does not capture the influence
of topological asymmetries. To tackle this issue, we analyze in more detail
the topological features of networks exhibiting solitary states in the
subsequent sections.

Figure 3: Map of regimes in the ($\sigma_{u},\sigma_{v}$)-parameter plane.
Color-coded is the probability of obtaining a solitary state for the following
topologies: (a) A nonlocally coupled ring network with $R=9$; (b) A random
(Erdős-Rényi) network with average degree $p=0.18$; (c) A scale-free
(Barabási-Albert) network with $m=10$. The probability is estimated from $100$
network realizations with different ICs uniformly distributed over the
intervals $u_{i}\in[-2.2,2.2]$, $v_{i}\in[-1.1,1.1]$. The network size is
$N=100$ in all cases.

To complement our analysis of the coupling intensities, we point out that in a
recent work for a globally coupled system (symmetric), Schülen et al. have
obtained similar patterns in the ($\sigma_{u},\sigma_{v}$)-parameter plane
[20]. For a class of such systems, the authors have derived the bifurcation
curves delimitating the regions of high probability for the occurrence of
solitary states [20]. The similarity of the patterns observed for the globally
coupled system and the networks studied here suggests that the same
bifurcation scenario gives rise to solitary states in the networks with
complex topologies. This fact is also investigated subsequently in this work.

##  3 The role of the network connectivity and local topological features

The results conveyed in Sec. 2 show that solitary states arise in networks
with complex topologies for the same range of coupling strength as for simple
regular networks. This raises a general question: What would be the
topological requirements for the onset of the solitary states? In order to
investigate the existence of possible constraints, we concentrate our efforts
on networks with topologies represented by BA graphs as they constitute the
example with asymmetric distribution of node degrees. As discussed in Sec. 2,
the connectivity of BA graphs is specified by the number of links $m$ added at
every step of the growth algorithm. With this, we investigate the onset of
solitary states in such networks for different levels of connectivity by
varying $m$. For that, we first adopt a statistical approach by considering
$1000$ realizations of the system Eqs. (1) with different ICs chosen randomly
as discussed in Fig. 3. In the top panel of Fig. 4(a), we show the share of
ICs resulting in solitary states with different number of solitary nodes for
varying $m$ as follows: one (green), two (orange), three (blue), four (orchid)
and five or more (red). The absence of solitary nodes (regime of
synchronization) is marked in gray. The colored region is the relative
fraction that of all ICs that goes to this specific state. For example, if
$m=10$, about $35\%$ of the simulations go to a state with a single solitary
oscillator, about $60\%$ go to either a single or two solitary nodes, so the
relative share of two solitaries is about $25\%$. Thus, they are given by the
differences between two curves. In effect, we have a $100\%$ stacked bar for
each value of $m$. For low values of $m$ (low connectivity), we observe that
solitary states occur for a small share of the ICs and they contain only one
solitary node. By increasing $m$, the share of realizations corresponding to
the completely synchronized network (gray) diminishes giving room to a larger
variety of solitary states. For instance, $m\approx 8$, the solitary states
are prevalent among all the realizations being observed with a different
number of solitary nodes. Despite the fact that they are observed for all
levels of network connectivity, the findings in Fig. 4(a) indicate that high
connectivity favors the formation of such states.

Figure 4: (a) Top panel: Share of network realizations yielding a solitary
state with one (green), two (orange), three (blue), four (light purple), and
five or more (red) solitary nodes for $m\in[1,25]$. The gray color corresponds
to realizations with no solitary states. Bottom panel: Ratio
$\langle.\rangle_{sol}/\langle.\rangle_{all}$ between average network
measures, node degree (blue), centrality (red), and average neighbor degree
(yellow). Parameters are: $\sigma_{u}=0.12$, $\sigma_{v}=0.12$, and $N=100$.
(b) Conditional probability $p(s|d)$ to obtain a solitary state for a given
node degree $d$ in a Barabási-Albert network with $m=4$ (green) and $m=8$
(magenta). (c) Conditional probability $p(s|d_{av})$ to obtain a solitary
state for a given average neighbor degree $d_{av}$ in a Barabási-Albert
network with $m=4$ (green) and $m=8$ (magenta). Parameters in (b) and (c):
$\sigma_{u}=0.12$, $\sigma_{v}=0.15$, and $N=500$.

Next, we refine our knowledge by analyzing how higher connectivity influences
the formation of solitary states. Specifically, we investigate local node
properties to find out if the success of high network connectivity can be
captured by measures such as the node degree, average neighbor node degree,
and eigenvector centrality. As the node degree $d_{i}$ has been already
defined in Sec. 2, the average neighbor node degree of node $i$ is given by
$d_{av,i}=\frac{\sum_{j=1}^{N}A_{ij}d_{j}}{d_{i}}$. This measure accounts for
the average level of connectivity of adjacent nodes of $i$. The eigenvector
centrality of a node $i$ is given by
$c_{i}=\frac{1}{\lambda}\sum_{j=1}^{N}A_{ij}c_{j}$, where $\lambda$ is the
largest eigenvalue calculated via $\mathbf{A}\mathbf{c}=\lambda\mathbf{c}$.
This measure ranks the node $i$ following the connectivity of its neighbors.
With this, we first obtain the average of these three measures over the entire
network, denoted by $\langle.\rangle_{all}$. Next, we obtain the average only
over the solitary nodes, which we denote by $\langle.\rangle_{sol}$. We
calculate the ratio $\langle.\rangle_{sol}/\langle.\rangle_{all}$ as a
function of the parameter $m$ in the bottom panel of Fig. 4(a). Interestingly,
for $m\lesssim 8$, the average over solitary nodes differs from the network
average. In particular, we point out the fact that the average neighbor node
degree is much higher for the solitary nodes, i.e., $\langle
d_{av,i}\rangle_{sol}>\langle d_{av,i}\rangle_{all}$ (yellow curve in the
bottom panel of Fig. 4(a)). This measure will be justified later on by the
bifurcation analysis (Sec. 4). In addition, by comparison with the top panel
of Fig. 4(a), we conclude that the asymmetries between the solitary nodes and
the nodes in the synchronized cluster favor the appearance of solitary states
with only one solitary node (green region in the top panel of Fig. 4(a)). In
contrast, at $m\approx 8$, the ratio
$\langle.\rangle_{sol}/\langle.\rangle_{all}$ of all measures approaches one,
indicating the establishment of homogeneity among the solitary and the
synchronized nodes. This configuration holds for $m>8$, giving rise to the
regime dominated by solitary states with a large number of solitary nodes.

To further analyze the dependencies observed in Fig. 4(a), we now investigate
the probability of a node being a solitary given one of its characteristics
discussed above such as node degree, average neighbor node degree, eigenvector
centrality. For that, we employ the concept of conditional probability, given
by

$p(s|x)=\frac{p(s\cap x)}{p(x)},$ (3)

where $p(s|x)$ stands for the probability of a node being a solitary $s$ given
a topological property $x$, e.g., the node degree $d$ or the average neighbor
node degree $d_{av}$. To calculate this conditional probability, we consider
$1000$ realizations with different ICs of BA networks possessing $N=500$ nodes
(for better statistical quality). With this, in the top panel of Fig. 4(b), we
obtain the conditional probability of a node being a solitary given its node
degree, i.e., $p(s|d)$. For $m=4$ (green circles), we observe that only nodes
with low $d<10$ degree have a nonzero probability of being solitary (see the
inset in the top panel of Fig. 4(b)). The probability of occurrence of
solitary nodes with a degree above this threshold is zero. This suggests that
the basins of attraction of the solitary states are rather small and does not
necessarily mean that solitary states are impossible for these node degrees.
On the other hand, for $m=8$ (pink circles), we find that both nodes at the
lower and the higher end of the degree distribution have a high probability of
becoming a solitary node. This inconsistency between the cases $m=4$ and $m=8$
suggests that the probability of a node being solitary is not only determined
by its respective degree. Therefore, we now look into the conditional
probability considering the average neighbor degree as the underlying
condition $p(s|d_{av})$. In Fig. 4(c) for both cases ($m=4$ and $m=8$), we
observe that a high average neighbor degree corresponds to a higher
probability of the respective node being solitary. This observation is in line
with the findings in the bottom panel of Fig. 4(a). Note that for $m\lesssim
8$ the average neighbor degree is the only measure that is on average higher
for the solitary nodes with respect to mean of the network, i.e. $\langle
d_{av,i}\rangle_{sol}>\langle d_{av,i}\rangle_{all}$. As $m$ increases, the
high connectivity in the network provides a higher average neighbor degree to
more nodes, creating the abundance of solitary states observed in the top
panel of Fig. 4(a).

##  4 Bifurcation analysis

The statistical analysis performed in Sec. 3 has shown that the average
neighbor node degree plays an important role in the formation of solitary
states. Now we deepen this knowledge by investigating the bifurcation scenario
giving rise to the solitary states in the considered complex network. To this
end, we make use of the fact that the entire synchronized cluster influences
the solitary node as a single input. That is, if node $s$ is a solitary and it
is connected to $d_{s}$ oscillators that are all in sync with each other, the
resulting coupling term in $s$ is given
$C_{s}=\sigma_{u}(u_{b}-u_{s})+\sigma_{v}(v_{b}-v_{s})$, where $(u_{b},v_{b})$
are the state variables of oscillators in the synchronized cluster. In turn,
an oscillator $n$, with node degree $d_{n}$, belonging to the synchronized
cluster, and directly coupled to the solitary node, has nonzero coupling input
only from the solitary node. The resulting coupling term in $n$ reads
$C_{n}=1/d_{n}\cdot[\sigma_{u}(u_{s}-u_{n})+\sigma_{v}(v_{s}-v_{n})]$. With
this, one can immediately infer that for large values of $d_{n}$, the coupling
term for the node $n$ vanishes and it decouples from the solitary node $s$.
Therefore, in order to ensure the coupling between the nodes $s$ and $n$, we
assume that $n$ has the minimum node degree $d_{n}$ among all neighbors of
$s$. In addition, we also assume that the rest of the network stays
synchronized at all times. With these assumptions, we can reduce the dynamics
of the network to the following equation of two interacting nodes [32]:

$\displaystyle\varepsilon\frac{du_{s}}{dt}$
$\displaystyle=f(u_{s},v_{s})+\sigma_{u}(u_{n}-u_{s})+\sigma_{v}(v_{n}-v_{s}),$
$\displaystyle\frac{dv_{s}}{dt}$ $\displaystyle=g(u_{s},v_{s}),$
$\displaystyle\varepsilon\frac{du_{n}}{dt}$
$\displaystyle=f(u_{n},v_{n})+\frac{1}{d_{n}}[\sigma_{u}(u_{s}-u_{n})+\sigma_{v}(v_{s}-v_{n})],$
$\displaystyle\frac{dv_{n}}{dt}$ $\displaystyle=g(u_{n},v_{n}),$ (4)

where the functions $f(u,v)$ and $g(u,v)$ are given in Eqs. (2). Following
this dimension reduction, we are now able to study the stability of solutions
of the system Eqs. (4). Considering the node degree $d_{n}$ as a continuous
parameter, we employ numerical continuation analysis to follow a solitary
solution as $d_{n}$ is varied. For this task, we use the software auto-07p
[33]. The resulting bifurcation diagram is shown in the top panel of Fig.
5(a). Even though, in reality the minimum node degree $d_{n}$ is an integer
number, the bifurcation analysis in Fig. 5(a) treating it as continuous
provides deep insight into the stability profile of the solitary node $s$.
Note that the limit-cycle hosting the trajectory of the solitary node appears
in a fold bifurcation at $d_{n}\approx 5.701$, see the red cross in the top
panel of Fig. 5(a). This analysis indicates that the neighbor nodes of a
solitary node must have a minimum node degree of $d_{n}\geq 6$. Therefore, we
establish this value of $d_{n}$ as a threshold value $d_{c}$ for the
occurrence of solitary states in the network. In addition, it is also visible
in the top panel of Fig. 5(a) that the increase of the parameter $d_{n}$ leads
to a period-doubling bifurcation transforming the dynamics of the solitary
node, see a blue cross. This observation suggests that the larger the neighbor
node degree, the higher is the period of the limit-cycle hosting the solitary
node.

Figure 5: (a) Bifurcation diagram of the reduced systems varying $d_{n}$ as a
bifurcation parameter. Red crosses indicate saddle-node bifurcations. The
vertical dashed line in magenta indicates the threshold value $d_{c}=7$. Top
panel: Reduced system in Eqs. (4). Bottom panel: Reduced system in Eqs. (5)
for $d_{s}=2$ (yellow), $d_{s}=4$ (red) and $d_{s}=8$ (purple). (b) The
minimal neighbor degree of all nodes in Barabási-Albert network with $m=4$
(top panel) and $m=8$ (bottom). Red circles indicate nodes that are observed
as solitaries. Blue circles are the nodes that are always in sync. The size of
the red circles is proportional to how often the corresponding node is
observed as a solitary. (c) Conditional probability $p(s|d_{min})$ to observe
a solitary for a given minimal neighbor node degree for Barabási-Albert
networks with $m=4$ (green) and $m=8$ (pink). Parameters: $\sigma_{u}=0.12$,
$\sigma_{v}=0.15$, and $N=500$.

The results shown in the top panel of Fig. 5(a) are obtained by assuming the
dynamics of two coupled oscillators as a descriptor of two interacting subsets
of the network, i.e., the synchronized cluster and the solitary node. This
approach can be improved by splitting the synchronized cluster into the
dynamics of a mean-field component and the dynamics of the node directly
coupled to the solitary node. With this, the dynamics of the solitary node $s$
and the directly coupled neighbor $n$ described in Eqs. (4) receive an extra
dynamical input corresponding to the synchronized mean-field. The dynamics of
the mean-field is described by an extra oscillator with state variables
$(u_{b},v_{b})$, which are not affected by the nodes $n$ and $s$. The
resulting equations for this version of the reduced system are given by:

$\displaystyle\varepsilon\frac{du_{b}}{dt}$ $\displaystyle=$ $\displaystyle
f(u_{b},v_{b}),\qquad\frac{dv_{b}}{dt}=g(u_{b},v_{b}),$
$\displaystyle\varepsilon\frac{du_{n}}{dt}$ $\displaystyle=$ $\displaystyle
f(u_{n},v_{n})+(1-\frac{1}{d_{n}})[\sigma_{u}(u_{b}-u_{n})+\sigma_{v}(v_{b}-v_{n})]+\frac{1}{d_{n}}[\sigma_{u}(u_{s}-u_{n})+\sigma_{v}(v_{s}-v_{n})],$
$\displaystyle\frac{dv_{n}}{dt}$ $\displaystyle=$ $\displaystyle
g(u_{n},v_{n}),$ $\displaystyle\varepsilon\frac{du_{s}}{dt}$ $\displaystyle=$
$\displaystyle
f(u_{s},v_{s})+(1-\frac{1}{d_{s}})[\sigma_{u}(u_{b}-u_{s})+\sigma_{v}(v_{b}-v_{s})]+\frac{1}{d_{s}}[\sigma_{u}(u_{n}-u_{s})+\sigma_{v}(v_{n}-v_{s})],$
$\displaystyle\frac{dv_{s}}{dt}$ $\displaystyle=$ $\displaystyle
g(u_{s},v_{s}),$ (5)

with $f(u,v)$ and $g(u,v)$ as in Eqs. (2). As before, the bifurcation analysis
is performed by varying the minimum degree $d_{n}$ of the neighbor node, see
the bottom panel of Fig. 5(a). We repeat the analysis for three different
values of the parameter attributed to the degree of the solitary node, namely
$d_{s}\in[2,4,8]$. Surprisingly, the choice of $d_{s}$ does not influence much
the threshold for the existence of solitary states. In all three cases, we
find the fold bifurcation located between $6<d_{n}<7$, slightly higher
compared to the case of the two-node reduction. This implies that the lower
threshold is $d_{n}\geq 7$. We can therefore confirm the claim that the
crucial parameter for the formation of solitary states is the minimal neighbor
node degree $d_{n}$. Neither the node degree of the solitary node nor the
_average_ neighbor node degree by itself have a causal role in the formation.

Finally, we statistically verify this result in the full system, i.e., via an
ensemble of realizations of the complex network. For that, we consider
realizations of BA networks with two different values of the parameter $m$ to
compute the quantity $d_{min,i}$, the minimum neighbor node degree of each
network node $i$. For $m=4$, in the top panel of Fig. 5$(b)$, we obtain
$d_{min,i}$ for all nodes $i$. In this figure, among the different
realizations, the nodes that are solitary at least once are colored in red,
while the other nodes are marked in blue. In addition, the size of the red
circles is proportional to how often the corresponding node is observed as a
solitary. We observe that the nodes with a higher value of $d_{min,i}$ are
more often solitaries, indicating that indeed the higher value of the minimum
neighbor node degree yields a higher probability of being solitary. Moreover,
for the $m=4$, the BA network has low connectivity implying that most nodes
have a $d_{min}$ below the threshold $d_{c}=7$ unraveled by the bifurcation
analysis in Fig. 5(a). The dashed horizontal line in Fig. 5(b) marks this
threshold and indeed all nodes below this threshold do not become solitary
(blue circles). Conversely, in a BA network with $m=8$, all nodes have a
degree larger than the critical threshold $d_{c}$. In such a network, all
nodes meet the criteria we established via bifurcation analysis. The network
simulations shown in the bottom panel of Fig. 5(b) confirm this hypothesis.
Note that $d_{min}$ for all nodes falls above the dashed horizontal line
marking $d_{c}$ and most of the nodes can be solitaries (red circles). Also,
in this case, higher $d_{min}$ corresponds to a higher probability of being
solitary. Finally, in Fig. 5(c), we complement our analysis by estimating the
conditional probability $p(s|d_{min})$ of a node $i$ with a given $d_{min,i}$
being solitary. For $m=4$ (green circles), the observed probability is nonzero
only for the nodes having $d_{min,i}>d_{c}$. In addition, the probability
increases for nodes with higher values of $d_{min,i}$. For $m=8$ (pink
circles), as all nodes already fall into the criteria for solitary
($d_{min,i}>d_{c}$ $\forall$ $i\in[1,500]$), there is no node with zero
probability of being a solitary. Again, the probability is higher for nodes
with higher $d_{min,i}$. Both cases confirm the existence of the threshold
$d_{c}$ for the minimum node degree of adjacent nodes provided by the
bifurcation analysis.

##  5 Discussion

In summary, by approaching the onset of solitary states in complex networks,
we have demonstrated that the connectivity of nodes adjacent to a solitary one
is an essential topological feature for the appearance of solitary states.
More specifically, we found a threshold for the minimum value of the degree of
nodes neighboring the solitary ones in the network. This finding is
statistically demonstrated by estimating the conditional probability of its
occurrence in ensembles of realizations of scale-free networks. Furthermore, a
dimensionality reduction of the network dynamics made possible a bifurcation
analysis confirming the existence of the minimum threshold.

Due to the symmetry-broken character of solitary states, their occurrence is
undesirable in many contexts. Therefore, the knowledge of topological
constraints for the onset of these states sheds light on a perspective of
their control, i.e., suppression or initiation. Finally, we emphasize that the
results reported here can be extended to other classes of networked systems
since there are no particular dynamical, or topological, restrictions on them.

Code and data are accessible upon request.

We thank Matthias Wolfrum and Alexander Gerdes for fruitful discussions. E.S.M
acknowledges the support by the Deutsche Forschungsgemeinschaft (DFG) via the
project number FE-359/22-1. This work was supported by the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) - Projektnummer -
163436311 - SFB 910.

## References

  * [1] S. Boccaletti, V. Latora, Y. Moreno, M. Chavez, and D.-U. Hwang.  Complex networks: Structure and dynamics.  Phys. Rep., 424(4-5):175–308, 2006.
  * [2] L. M. Pecora, F. Sorrentino, A. M. Hagerstrom, T. E. Murphy, and R. Roy.  Cluster synchronization and isolated desynchronization in complex networks with symmetries.  Nat. Comm., 5(1):1–8, 2014.
  * [3] W. Poel, A. Zakharova, and E. Schöll.  Partial synchronization and partial amplitude death in mesoscale network motifs.  Phys. Rev. E, 91(2):022915, 2015.
  * [4] D. M. Abrams and S. H. Strogatz.  Chimera states for coupled oscillators.  Phys. Rev. Lett., 93:174102, Oct 2004.
  * [5] A. Zakharova.  Chimera Patterns in Networks: Interplay between Dynamics, Structure, Noise and Delay.  Springer, 2020.
  * [6] Y. Maistrenko, B. Penkovsky, and Rosenblum M.  Solitary state at the edge of synchrony in ensembles with attractive and repulsive interactions.  Phys. Rev. E, 89:060901, 2014.
  * [7] P. Jaros, Y. Maistrenko, and T. Kapitaniak.  Chimera states on the route from coherence to rotating waves.  Phys. Rev. E, 91(2):022907, 2015.
  * [8] P. Jaros, S. Brezetsky, R. Levchenko, D. Dudkowski, T. Kapitaniak, and Y. Maistrenko.  Solitary states for coupled oscillators with inertia.  Chaos, 28(1):011103, 2018.
  * [9] F. Hellmann, P. Schultz, P. Jaros, R. Levchenko, T. Kapitaniak, J. Kurths, and Y. Maistrenko.  Network-induced multistability through lossy coupling and exotic solitary states.  Nat. Comm., 11(592), 2020.
  * [10] E. Rybalova, N. Semenova, G.I. Strelkova, and V.S. Anishchenko.  Transition from complete synchronization to spatio-temporal chaos in coupled chaotic systems with nonhyperbolic and hyperbolic attractors.  Eur. Phys. J. Special Topics, 226(9):1857–1866, 2017.
  * [11] E.V. Rybalova, G.I. Strelkova, and V.S. Anishchenko.  Mechanism of realizing a solitary state chimera in a ring of nonlocally coupled chaotic maps.  Chaos, Solitons & Fractals, 115:300–305, 2018.
  * [12] N. Semenova, T. Vadivasova, and V.S. Anishchenko.  Mechanism of solitary state appearance in an ensemble of nonlocally coupled lozi maps.  Eur. Phys. J. Special Topics, 227(10):1173–1183, 2018.
  * [13] I. Franović, S. Eydam, N. Semenova, and A. Zakharova.  Unbalanced clustering and solitary states in coupled excitable systems.  Chaos, 32(1):011104, 2022.
  * [14] M. Mikhaylenko, L. Ramlow, S. Jalan, and A. Zakharova.  Weak multiplexing in neural networks: Switching between chimera and solitary states.  Chaos, 29(2):023122, 2019.
  * [15] L. Schülen, D.A. Janzen, E.S. Medeiros, and A. Zakharova.  Solitary states in multiplex neural networks: Onset and vulnerability.  Chaos, Solitons & Fractals, 145:110670, 2021.
  * [16] S. Majhi, T. Kapitaniak, and D. Ghosh.  Solitary states in multiplex networks owing to competing interactions.  Chaos, 29(1):013108, 2019.
  * [17] E. Rybalova, A. Zakharova, and G.I. Strelkova.  Interplay between solitary states and chimeras in multiplex neural networks.  Chaos, Solitons & Fractals, 148:111011, 2021.
  * [18] L. Schülen, S. Ghosh, A.D. Kachhvah, A. Zakharova, and S. Jalan.  Delay engineered solitary states in complex networks.  Chaos, Solitons & Fractals, 128:290 – 296, 2019.
  * [19] E. Rybalova, V. S. Anishchenko, G. I. Strelkova, and A. Zakharova.  Solitary states and solitary state chimera in neural networks.  Chaos, 29(7):071106, 2019.
  * [20] L. Schülen, A. Gerdes, M. Wolfrum, and A. Zakharova.  The solitary route to chimera states.  arXiv preprint arXiv:2204.00385, 2022.
  * [21] T. Nishikawa and A. E. Motter.  Symmetric states requiring system asymmetry.  Phys. Rev. Lett., 117:114101, Sep 2016.
  * [22] J. D. Hart, Y. Zhang, R. Roy, and A. E. Motter.  Topological control of synchronization patterns: Trading symmetry for stability.  Phys. Rev. Lett., 122:058301, 2019.
  * [23] Everton S. Medeiros, Ulrike Feudel, and Anna Zakharova.  Asymmetry-induced order in multilayer networks.  Phys. Rev. E, 104:024302, 2021.
  * [24] J. F. Yang, T. A. Berrueta, A. M. Brooks, A. T. Liu, G. Zhang, D. Gonzalez-Medrano, S. Yang, V. B Koman, P. Chvykov, M. Z. Miskin, et al.  Emergent microrobotic oscillators via asymmetry-induced order.  arXiv preprint arXiv:2205.09814, 2022.
  * [25] I. Omelchenko, O. Omel’chenko, P. Hövel, and E. Schöll.  When nonlocal coupling between oscillators becomes stronger: Patched synchrony or multi-chimera states.  Phys. Rev. Lett., 110:224101, 05 2013.
  * [26] D.W. Storti and P.G. Reinhall.  Phase-locked mode stability for coupled van der pol oscillators.  J. Vib. Acoust., 122(3):318–323, 2000.
  * [27] L.A. Low, P.G. Reinhall, and D.W. Storti.  An investigation of coupled van der pol oscillators.  J. Vib. Acoust., 125(2):162–169, 2003.
  * [28] L.A. Low, P.G. Reinhall, D.W. Storti, and E.B. Goldman.  Coupled van der pol oscillators as a simplified model for generation of neural patterns for jellyfish locomotion.  Struct. Control Health Monit., 13(1):417–429, 2006.
  * [29] B. Bollobás.  Random graphs.  In Modern graph theory, pages 215–252. Springer, 1998.
  * [30] R. Albert and A.-L. Barabási.  Statistical mechanics of complex networks.  Rev. Mod. Phys., 74(1):47, 2002.
  * [31] A.-L. Barabási.  Network science.  Philos. Trans. R. Soc., 371(1987):20120375, 2013.
  * [32] M. Golubitsky and I. Stewart.  The symmetry perspective: from equilibrium to chaos in phase space and physical space, volume 200.  Springer Science & Business Media, 2003.
  * [33] E. J. Doedel, A. R. Champneys, F. Dercole, T. F. Fairgrieve, Y. A. Kuznetsov, B. Oldeman, R.C. Paffenroth, B. Sandstede, X.J. Wang, and C.H. Zhang.  Auto-07p: Continuation and bifurcation software for ordinary differential equations, 2007.


# Neural network as a tool for design of amorphous metal alloys with desired
elastoplastic properties

B.N. Galimzyanov<EMAIL_ADDRESS>M.A. Doronina
<EMAIL_ADDRESS>A.V. Mokshin<EMAIL_ADDRESS>Kazan Federal
University, 420008 Kazan, Russia

###### Abstract

The development and implementation of the methods for designing amorphous
metal alloys with desired mechanical properties is one of the most promising
areas of modern materials science. Here, the machine learning methods appear
to be a suitable complement to empirical methods related to the synthesis and
testing of amorphous alloys of various compositions. In the present work, it
is proposed a method to determine amorphous metal alloys with mechanical
properties closest to those required. More than $50\,000$ amorphous alloys of
different compositions have been considered, and the Young’s modulus $E$ and
the yield strength $\sigma_{y}$ have been evaluated for them by the machine
learning model trained on the fundamental physical properties of the chemical
elements. Statistical treatment of the obtained results reveals that the
fundamental physical properties of the chemical element with the largest mass
fraction are the most significant factors, whose values correlate with the
values of the mechanical properties of the alloys, in which this element is
involved. It is shown that the values of the Young’s modulus $E$ and the yield
strength $\sigma_{y}$ are higher for amorphous alloys based on Cr, Fe, Co, Ni,
Nb, Mo and W formed by the addition of semimetals (e.g. Be, B, Al, Sn),
nonmetals (e.g. Si and P) and lanthanides (e.g. La and Gd) than for alloys of
other compositions. Increasing the number of components in alloy from $2$ to
$7$ and changing the mass fraction of chemical elements has no significantly
impact on the strength characteristics $E$ and $\sigma_{y}$. Amorphous metal
alloys with the most improved mechanical properties have been identified. In
particular, such extremely high-strength alloys include Cr80B20 (among
binary), Mo60B20W20 (among ternary) and Cr40B20Nb10Pd10Ta10Si10 (among
multicomponent).

###### keywords:

machine learning; materials design; mechanical properties; metals; amorphous
alloys

##  1 Introduction

Amorphous metal alloys are the promising materials for the automotive,
aerospace, energy, electronics and medical technology industries [1, 2, 3, 4].
High corrosion resistance, high magnetic permeability, superior mechanical
strength, high fracture toughness, high elastic strain limit and high
formability are just some of the unique set of properties that make amorphous
metal alloys widely applicable [5, 6, 7]. Such the combination of properties
is directly due to the absence of structural order accompanied by defects that
is typical for crystalline analogues [8, 9, 10, 11]. However, despite all the
advantages of amorphous metal alloys, their production is complicated by the
fact that the formation of a stable disordered structure depends strongly on
alloy composition (i.e. number of components, type of added chemical elements)
and its preparation protocol (i.e. cooling and compression procedures, initial
and final melt temperatures) [12, 13, 14, 15, 16].

Amorphous metal alloys are actively studied for more than 80 years, beginning,
in particular, with the works of Kramer [17, 18]. One of the first methods of
practical formation of alloys with amorphous structure was based on the so-
called electrodeposition process. Later, in the 60’s of the 20th century, the
first works related with formation of amorphous metal films by rapid cooling
of the corresponding melts were appeared [19, 20]. As it turned out later,
amorphization of metallic melts of almost any composition is possible if
extremely fast cooling is used. The next stage in the development of the
amorphous alloy formation methodology concerned the consideration of alloys in
eutectics, where it was found that bulk amorphous samples of more than $1$ mm
thickness can be formed [21]. Further attention in this area focuses on some
aspects. Namely, the mechanical properties of bulk amorphous metal alloys are
strongly dependent on alloy composition and chemical purity of raw material.
The strength properties of amorphous alloys can be significantly reduced due
to the presence of impurities. Moreover, bulk amorphous metal alloys are
inherently fragile. Therefore, in the early 2000’s, studies were aimed to
improve the alloy hardening methods as well as to determine the relationship
between the key mechanical properties of amorphous metal alloys, which include
the Young’s modulus $E$, the yield strength $\sigma_{y}$ and the strength
$\sigma_{f}$ [22, 23]. It has been shown that the relationship between the
hardness $H$ (by Vickers method), the strength $\sigma_{f}$, the Young’s
modulus $E$ and the yield strength $\sigma_{y}$ of amorphous metal alloys is
close to linear and can be reproduced, for example, by Tabor’s relation
$H=K\sigma_{y}$, by Johnson’s model $H=\sigma_{y}(a+b\ln[cE/\sigma_{y}])$ and
by relation $\sigma_{f}=dE^{1/2}$ (here, $K$, $a$, $b$, $c$ and $d$ are
constants) [24, 25, 26]. These studies found that amorphous metal alloys with
large values of $E$ and $\sigma_{y}$ are characterized by high hardness $H$
and strength $\sigma_{f}$.

The synthesis of amorphous metal alloy with desired mechanical properties may
require listing various combinations of compositions followed by mechanical
testing. This makes the process of synthesizing new alloys extremely difficult
and significantly increases the costs. Then, methods of computer design seem
to be a suitable support for empirical methods at the stage of determining
amorphous metal alloys with desired mechanical properties [27, 28]. In recent
decades, rapid development of information technologies as well as automation
of data collection and storage processes contribute to accumulation and
systematization of information about the physical and mechanical properties of
bulk amorphous metal alloys glasses [29, 30, 31, 32]. The methods of machine
learning operate with large arrays of the data and allow us to determine the
relationship between composition and properties of alloys both already known
and not previously known [33, 34, 35, 36]. For example, Xiong and co-authors
have been developed a machine learning model that can predict the glass-
forming ability and elastic moduli of bulk metallic glasses based on the
fundamental atomic properties, chemical and physical properties obtained from
experiments or density functional theory simulations [37]. These results find
the importance of individual chemical element properties and macroscopic
properties in determining the strength characteristics of amorphous alloys.
The results obtained by Khakurel et al. established that the average
concentration of valence electrons, the atomic radius and the melting
temperature are the key properties, which are correlated with the Young’s
modulus of compositionally complex alloys [38]. The results of this work can
also be extended to amorphous metal alloys, as it is confirmed in Refs. [39,
40]. In addition, as it was found in Ref. [41] using a machine learning model,
the Young’s modulus of metal alloys under normal conditions correlates with
the yield strength and with the glass transition temperature. In this case,
the specificity of “chemical formula” of alloy, which is determined by the
molar mass and the number of components, is not as important as is usually
expected. Johnson and Samwer have found that the mechanical properties
(elastic constants, compressive yield strength, elastic strain limit) of $30$
bulk metallic glasses as functions of the scaled temperature $T_{R}/T_{g}$
obey the universal law $\propto a-b(T_{R}/T_{g})^{2/3}$, where $a$ and $b$ are
the constants, $T_{R}$ is the room temperature, $T_{g}$ is the glass
transition temperature [42]. The results of this work systematize existing
knowledge about the mechanical properties of amorphous alloys. An artificial
neural network has created by Jeon and co-authors for designing Fe-based
amorphous metal alloys with the desired crystallization temperature and glass
transition temperature [43]. Thus, all these studies show that the machine
learning methods are suitable tool to find new amorphous alloys with required
physical and mechanical properties. Despite the significant number of such
studies, little attention has been paid to the development of methods for
determining previously unknown amorphous alloys with the desired mechanical
properties.

The present work proposes a new method for determining amorphous metal alloys
of arbitrary composition based on a large set of empirical data. The
originality of this method is that it is based on a machine learning model
capable of predicting the Young’s modulus and the yield strength of amorphous
alloys taking into account the fundamental properties of each chemical element
that forms the alloys. It is quite significant that the obtained results lead
to new knowledge, which will contribute to the determination of amorphous
metal alloys that maximally satisfy the required mechanical properties.

##  2 Method for determining the mechanical properties of amorphous metal
alloys

###  2.1 General strategy of the method

The developed method for determining amorphous metal alloys is based on a
machine learning model, which is an artificial neural network of direct
propagation. The main advantage of this method is the possibility to calculate
the Young’s modulus $E$ and the yield strength $\sigma_{y}$ both for known
amorphous metal alloys and for alloys that are yet to be synthesized. The
developed method makes it possible to determine $E$ and $\sigma_{y}$ of
alloys, whose number of components varies in the range from $2$ to $7$. Note
that such the number of components is ordinary for the majority of known metal
alloys. In addition, the proposed method can be adapted to identify alloys
with large number of components at the presence of appropriate data for neural
network training. The composition and mass fraction of chemical elements in
the generated alloys are the control parameters, which allow us to construct a
diverse set of compounds.

Figure 1: Four-stage scheme of the method for determining amorphous metal
alloys and calculating their mechanical properties.

The general strategy for determining amorphous metal alloys implemented in
this work consists of four main stages [see Figure 1]:

  * 1.

Stage I. This stage includes the process of data collection and
systematization of information about the properties of multicomponent
amorphous metal alloys based on Al, Au, Ca, Co, Cu, Fe, La, Hf, Mg, Ni, Pd,
Pt, Sc, Ti, W, Zr, etc., as well as information about the properties of the
other additional chemical elements involved in the formation of these alloys.
Among these properties are the atomic mass $m_{a}$, the covalent radius
$r_{c}$, the ionization energy $E_{i}$ and the electronegativity $\chi$, which
characterize the nature of the chemical element [see Table 1]. This choice is
due to the following reasons. First, these parameters most clearly define the
possible physical and chemical bonds between the elements, which can either
promote or inhibit the formation of an amorphous structure. For example,
according to the empirical rule proposed by Inone et al. in the early 1990’s
[44], the difference in atomic sizes must be greater than $12$ % for good
amorphization of a liquid. Secondly, most of the intrinsic properties of
chemical elements (especially of the same type) are correlated. In addition,
the thermal conductivity $\lambda$, the specific heat capacity $C_{s}$, the
density $\rho$, the melting temperature $T_{m}$ and the boiling temperature
$T_{b}$ of chemical elements at normal conditions are used. The atomic number
$Z$ and the mass fraction $m_{f}$ of each chemical element in the alloy are
used to characterize the alloy composition. The Young’s modulus $E$ and the
yield strength $\sigma_{y}$ are also applied, whose values are known for the
considered amorphous alloys. The values of all the listed physical properties
are taken from the database ITPhyMS (Information technologies in physical
materials science) [45] and the database Materials Project [46] as well as
from Refs. [36, 47, 48, 49, 50] [see Supplementary data of the present work].
These properties are characterized by different physical dimensions and by
different ranges of values. Therefore, the properties are calibrated so that
their values vary in the range [$0$; $1$]. The calibration is done according
to the rule

$\text{Property}^{\prime}=\frac{\text{Property}-\text{Value}_{\text{min}}}{\text{Value}_{\text{max}}-\text{Value}_{\text{min}}},$
(1)

where “$\text{Value}_{\text{min}}$” and “$\text{Value}_{\text{max}}$” are the
smallest and largest known values of the “Property”. Moreover, all these
listed properties correlate with the mechanical properties of materials. For
example, Xiong et al. have shown that the accuracy of predicting the
mechanical properties of amorphous metal alloys is improved when the
quantities $m_{a}$, $Z$, $r_{c}$, $\rho$, $\lambda$, $T_{m}$ and $T_{b}$ are
considered in a machine learning model [37]. In addition, the results obtained
by Wang based on the analysis of a large set of empirical data for amorphous
alloys allow one to establish the existence of correlation between elastic
moduli (i.e. Young’s modulus, shear modulus, bulk modulus), microstructural
features, rheological properties, the glass transition temperature, the
melting temperature and the boson peak [47, 48].

  * 2.

Stage II. Alloys with different compositions are generated. Taking into
account the number of possible components, combinations of all chemical
elements and their mass fraction, up to $10^{18}$ different compositions can
be determined simultaneously. When obtaining alloys, those chemical elements
are selected that are included on the alloys in the training dataset. In the
present work, $32$ chemical elements were used including transition metals
(Fe, Co, Ni, Cu, etc.), semimetals (B, Al, Sn, etc.), lanthanides (La, Gd, Er,
etc.) and alkali and alkaline earth metals (Li, Be, Mg, Ca, etc.). A list of
all the considered chemical elements is given on Table S1 in Supplementary
data. The mass fraction of the chemical elements in a generated alloy is also
set randomly so that the total mass fraction of all chemical elements is equal
to $100$ %. A set of physical properties is created for each chemical element
[see Table 1].

  * 3.

Stage III. Information about the alloy composition and the physical properties
of all the chemical elements is processed by the pre-trained neural network.
This neural network evaluates the Young’s modulus $E$ and the yield strength
$\sigma_{y}$ for all generated alloys. The training procedure of the neural
network is discussed in more details in the subsection “Machine learning
model: structure and training”.

  * 4.

Stage IV. Statistical interpretation of machine learning results is performed.

Thus, the proposed method makes it possible to perform a complete cycle of
alloy design and determine its mechanical properties: from obtaining the
correct alloy composition to calculating the correct values of $E$ and
$\sigma_{y}$.

Table 1: Physical properties of chemical elements used as input parameters in the artificial neural network. Property | Symbol | Unit
---|---|---
Atomic number | $Z$ | –
Mass fraction of elements | $m_{f}$ | %
Atomic mass | $m_{a}$ | a.e.m.
Covalent radius | $r_{c}$ | pm
Ionization energy | $E_{i}$ | eV
Electronegativity | $\chi$ | –
Thermal conductivity | $\lambda$ | W/(m$\cdot$K)
Specific heat capacity | $C_{s}$ | J/(g$\cdot$K)
Density | $\rho$ | g/cm3
Melting temperature | $T_{m}$ | K
Boiling temperature | $T_{b}$ | K

###  2.2 Machine learning model: structure and training

The machine learning model is the four-layer artificial neural network. The
first layer has $77$ input neurons for the values of $11$ physical properties
for all chemical elements of the obtained alloy ($7$ input neurons are
allocated to each property because the maximal number of components in the
alloy is also seven). If the number of components in the alloy is less than
$7$, then the remaining neurons are unused. The next two layers are hidden.
The first hidden layer consists of $80$ neurons, while the second hidden layer
has $10$ neurons. Note that the number of neurons in the hidden layers is
optional. The neural network produces close results with $80$ to $100$ neurons
in the first hidden layer and with $10$ to $80$ neurons in the second hidden
layer. The fourth layer consists of one neuron that determines the Young’s
modulus $E$ or the yield strength $\sigma_{y}$. It is important to note that
two separate independent neural networks with the same structure are used to
calculate the values of $E$ and $\sigma_{y}$.

Calculation of the values of all neurons is carried out by expression [51]:

$n_{i}^{(k)}=f\left(\sum_{j=1}^{N_{k-1}}w_{ij}^{(k-1)}n_{j}^{(k-1)}+b_{i}^{(k)}\right).$
(2)

Here, $n_{i}^{(k)}$ is the value of the $i$th neuron in the $k$th layer
($k=2,\,3,\,4$); $w_{ij}^{(k-1)}$ is the value of the $(k-1)$th layer weight
going from a neuron with index $j$ to a neuron with index $i$ from the $k$th
layer; $b_{i}^{(k)}$ is the bias weight acting on a neuron with index $i$;
$N_{k-1}$ is the number of neurons in the $(k-1)$th layer. The sigmoid
$f(x)=1/(1-\exp[-x])$ is applied as the activation function [52].

Figure 2: (a) Plot of the predicted Young’s modulus $E$ versus the empirical
$E$. (b) Plot of the predicted yield strength $\sigma_{y}$ versus the
empirical $\sigma_{y}$. Top insets: mean relative error as function of the
number of training cycles for $E$ and $\sigma_{y}$. Bottom insets: dependence
of the mean relative error and the training cycles on the training rate
$\gamma$, from which the optimal value of $\gamma$ (indicated by the red
arrows) was determined.

The neural network is trained using the backpropagation algorithm [53, 54].
The values of the weight coefficients are adjusted as follows:

$w_{ij}^{(k),\,new}=w_{ij}^{(k)}-\gamma\frac{\partial\xi}{\partial
w_{ij}^{(k)}},$ (3)

where $\xi$ is the squared error between the output neuron and the desired
value of the mechanical property; $\gamma$ is the training rate. In the
present work, the training rate is $\gamma=0.3$, which is optimal for the
created neural network. At the training rate $\gamma=0.3$, the machine
learning model gives the best result for $E$ and $\sigma_{y}$ with the lower
MRE at the relatively small number of training cycles [see insets on Figures
2(a) and 2(b)]. The original dataset is divided into training and validation
subsets in proportion 80:20. The training subset consists of amorphous metal
alloys based on Al, Au, Ca, Co, Cu, Fe, La, Hf, etc. with different
compositions, for which the values of $E$ and $\sigma_{y}$ are known [36, 47,
48, 49, 50]. The physical properties of the chemical elements of these alloys
are also used in the training procedure [see Supplementary data]. To verify
the correctness of the machine learning results, the validation subset is
applied, which includes amorphous alloys that were not included on the
training subset. The criterion to stop the training procedure is the minimal
error between the results of the output neuron and the required values from
the validation subset.

###  2.3 Validation of the machine learning model

Typically, RMSRE (Root Mean Squared Relative Error), RMSE (Root Mean Square
Error), RRMSE (Relative Root Mean Square Error), MSE (Mean Square Error), MAE
(Mean Absolute Error) or MRE (Mean Relative Error) are used as indicators for
measuring accuracy of results [55, 56, 57, 58, 59]. In the present work, it
was important to use an indicator that does not depend on units of physical
quantities. At the same time, this indicator must be easy to estimate.
Therefore, we chose the MRE, which is calculated by expression:

$\text{MRE}=\frac{1}{N}\sum_{i=1}^{N}\frac{|\mathcal{M}_{\text{ANN}}-\mathcal{M}_{\text{req}}|}{\mathcal{M}_{\text{req}}}\times
100\%.$ (4)

Here, $\mathcal{M}=\\{E$ or $\sigma_{y}\\}$ denotes the mechanical property;
$\mathcal{M}_{\text{ANN}}$ is the result of the neural network;
$\mathcal{M}_{\text{req}}$ is the required value of the mechanical property;
$N$ is the number of items in the validation subset. We find that the MRE is
$\sim 13$ % for the Young’s modulus $E$ and $\sim 11$ % for the yield strength
$\sigma_{y}$. This relatively low MRE indicates a correlation between the
predicted and empirical values of the mechanical properties that is also
confirmed by the results presented in Figures 2(a) and 2(b). Moreover, the
values of the MRE are stable. This is confirmed by the computed loss functions
[see insets on Figures 2(a) and 2(b)], which reach a plateau after $3\times
10^{3}$ training cycles. Thus, the results of the machine learning model are
reliable and predictable.

##  3 Properties importance scores

The analysis of the importance scores shows that all the considered physical
properties ($\lambda$, $T_{b}$, $\chi$, $\rho$, $C_{s}$, $T_{m}$, $E_{i}$,
$r_{c}$ and $m_{a}$) are necessary for the correct evaluation of the Young’s
modulus $E$ and the yield strength $\sigma_{y}$ by the machine learning model.
As can be seen from Figure 3(a), in the case of the Young’s modulus $E$, the
importance scores of the properties $\lambda$, $T_{b}$, $C_{s}$, $T_{m}$ and
$m_{a}$ are similar (MRE$\sim 45$ %), and these physical properties can be
recognized as significant factors. The lowest importance scores are observed
for the density $\rho$ and for the properties characterizing the chemical
nature of the atoms: the ionization energy $E_{i}$ (MRE$\sim 75$ %), the
electronegativity $\chi$ (MRE$\sim 61$ %) and the covalent radius $r_{c}$
(MRE$\sim 60$ %). These properties have less impact on the result of the
machine learning model. Furthermore, considering all the properties reduces
the error to MRE$\sim 13$ %. The importance scores of the properties in the
case of the yield strength $\sigma_{y}$ differ significantly from the Young’s
modulus $E$ [see Figure 3(b)]. The main factors affecting the yield strength
$\sigma_{y}$ are the thermal conductivity $\lambda$ (MRE$\sim 22$ %) and the
covalent radius $r_{c}$ (MRE$\sim 23$ %). MRE for other parameters is above
$28$ %. Together, these properties lead to the best result, where the error is
MRE$\sim 11$ %.

Figure 3: Importance scores for each physical property: (a) in the case of the
Young’s modulus and (b) in the case of the yield strength. (c) Mean relative
error as a function of the number of parameters in the artificial neural
network input. (d) Pearson correlation heat map for the considered physical
properties.

Thus, the formation of the machine learning models for the Young’s modulus $E$
and the yield strength $\sigma_{y}$ using a single physical property
($\lambda$, or $T_{b}$, or $\chi$, or $\rho$, or $C_{s}$, or $T_{m}$, $\dots$)
produces an error that is much larger than the error when these machine
learning models are formed with the entire set of physical properties.
Mathematically, such the situation is possible when the correlation between
the parameter $E$ (or $\sigma_{y}$) and some individual physical property
appears indirectly (not explicitly). In turn, this means that it is not
possible to obtain analytical expressions relating a mechanical property with
any parameter of the set ($\lambda$, $T_{b}$, $\chi$, $\rho$, $C_{s}$,
$T_{m}$, $E_{i}$, $r_{c}$ and $m_{a}$) and correctly reproducing the results
for an arbitrary metal alloy. The methodology of artificial neural networks
used in this study makes it possible to obtain a correspondence between the
Young’s modulus $E$ (or the yield strength $\sigma_{y}$) and the whole set of
the considered physical properties; and this correspondence is reproduced not
by an analytical expression, but by the internal structure of the formed
neural network. In fact, this feature of this methodology is an advantage when
dealing with a fairly large set of parameters.

An additional evaluation of the accuracy of the machine learning model was
performed by computing the MRE for different numbers of physical properties in
the input of the neural network. Figure 3(c) shows that at adding properties
in the order of $\rho$, $T_{m}$, $T_{b}$, $\lambda$, $C_{s}$, $m_{a}$,
$r_{c}$, $\chi$ and $E_{i}$ the MRE decreases from $\sim 58$ % to $\sim 13$ %
for $E$ and from $\sim 33$ % to $\sim 11$ % for $\sigma_{y}$. A rapid decrease
of the error is observed when the temperatures $T_{m}$ and $T_{b}$ as well as
the quantities $E_{i}$ and $\chi$ have been added, which may be due to their
multicollinearity. Figure 3(d) shows that the Pearson correlation coefficients
for the considered properties take both positive and negative values in the
range from $-1$ to $1$ [60]. For example, the positive correlation between the
temperatures $T_{m}$ and $T_{b}$ is due to the fact that an increase in the
melting temperature leads to an increase in the boiling temperature [61]. An
increase in the atomic mass $m_{a}$ of the alloy components usually leads to
an increase in its density $\rho$, which leads to a positive correlation
between $m_{a}$ and $\rho$ [62, 63]. The presence of the pronounced negative
correlation between the pairs $r_{c}$, $E_{i}$ and $r_{c}$, $\chi$ is due to
the fact that a decrease in the covalent radius $r_{c}$ leads to an increase
$E_{i}$ and $\chi$ by increasing the electron density in the atom [64, 65].

##  4 Statistical interpretation of the results

In the present study, $50\,000$ different amorphous metal alloys were obtained
by the proposed method. All alloys were sorted according to the atomic number
$Z$ of the basic chemical element and the number of components in the alloy.
Using the trained machine learning model, the values of $E$ and $\sigma_{y}$
were calculated for each alloy. Then, the average value of the mechanical
property was found for each $X$-based alloy consisting $n$ components (where
$n=2,\,3,...,\,7$). Here, $X$ denotes a chemical element, the mass fraction of
which in the alloy is greater than that of other elements. For example, Al-
based binary alloys Al90Fe10, Al80Cu20, Al60Ni40, etc. were selected and the
average values of $E$ and $\sigma_{y}$ were determined for all these alloys.
Similar calculations were performed for alloys based on other metals with
different number of components. Then, the dependence of the average values of
$E$ and $\sigma_{y}$ on the atomic number $Z$ of the basic chemical element
has been determined.

Figure 4: 3D plot of the dependence of the mechanical properties on the atomic
number $Z$ of the basic chemical element in the alloy and on its number of
components: (a) for the Young’s modulus $E$ and (b) for the yield strength
$\sigma_{y}$.

In the statistical interpretation, the results reveal that $E$ and
$\sigma_{y}$ depend mainly on the properties of the chemical element with the
largest mass fraction. As seen in Figure 4, changing the number of components
in alloy has no significant effect on values of $E$ and $\sigma_{y}$. In the
array of $50\,000$ different alloys obtained by the machine learning model,
some alloys with the highest Young’s modulus $E$ and the yield strength
$\sigma_{y}$ were selected [see Table 2]. The results show that these alloys
are mainly based on Ti, Cr, Fe, Co, Ni, Zr, Nb, Mo, Pd, Ta and W. For example,
these are Mo60B20W20, Co40B20Be20Al20, Cr40B20Nb10Pd10Ta10Si10 and
Cr30Mo20W20Pd10Gd10B10 alloys for which the mechanical properties are $E>300$
GPa and $\sigma_{y}>5.0$ GPa. It is important to note that the results for
alloys of these compositions were not previously known, although alloys of
some related compositions have been studied. So, for example, for W46Ru37B17,
Co43B31.5Fe20Ta5.5 and Co60B35Ta5 it was experimentally established that
$E>250$ GPa and $\sigma_{y}>5.0$ GPa [50, 41]. Obtained results reveal that
the alloys based on Cr, Mo and W from the group VI-B of the Periodic Table of
the Elements have improved mechanical properties. Such the metals as Cr, Mo
and W are refractory and have very high hardness [66, 67, 68]. Then, their
significant presence in an alloy improves its strength. Note that this fact is
also known in metallurgy, where these metals are widely used to increase the
hardness of steel alloys, to increase wear resistance and to form wear-
resistant coatings (e.g. alloys Cr-Co, Cr-Fe, Mo-Fe, Mo-Cr-Fe, W-Fe, W-Ni-Co)
[69, 70]. The machine learning model predicts improved mechanical properties
in the case of alloys based on Ti, Cr, Fe, Co, Ni, Zr, Nb, Mo, etc. when these
alloys are doped with other metals (e.g. Be, B, Hf), nonmetals (e.g. Si, P)
and lanthanides (e.g. La, Gd). The mechanical properties of alloys based on
Al, Mg, Ca, Cu, Zn, Ag, Au, Hf, lanthanides, etc. are inferior to those of the
alloys based on Ti, Cr, Fe, Co, etc. It should be noted that the relatively
high values of $E$ and $\sigma_{y}$ purely for B and Si are due to statistical
error, since B-based and Si-based alloys were not used in the training stage
of the machine learning model. However, alloys containing B and Si were
considered.

Table 2: Young’s modulus $E$ and yield strength $\sigma_{y}$ predicted by machine learning model for different amorphous metal alloys. Here, the average accuracy is $\sim 88$ %. Number of components | Alloy |  $E$, GPa | Alloy |  $\sigma_{y}$, GPa
---|---|---|---|---
2 |  Cr80B20 | $305$ |  Pd60B40 | $5.13$
|  W60Hf40 | $271$ |  W60Hf40 | $5.04$
3 |  Mo60B20W20 | $319$ |  Mo60B20Si20 | $5.35$
|  Nb40Hf40B20 | $289$ |  Ni60B20Sc20 | $5.17$
|  Ni60B20W20 | $280$ |  Pd60B20P20 | $4.31$
4 |  Ag40B20Sc20Ta20 | $302$ |  Co40B20Be20Al20 | $5.27$
|  Zr40Ni20B20Be20 | $285$ |  Nb40W20La20B20 | $4.88$
|  Cr40B20Zr20Hf20 | $271$ |  Ti40W20Pd20B20 | $4.70$
5 |  Ti30Fe30B20Sn10Be10 | $296$ |  Co40B30Ag10Gd10Si10 | $5.54$
|  Pd40B20Si20P10Hf10 | $289$ |  Fe50B20Mo10Ta10Ag10 | $5.39$
6 |  Cr40B20Nb10Pd10Ta10Si10 | $310$ |  Cr30Mo20W20Pd10Gd10B10 | $5.62$
|  Pd40Be20Mo10Ti10B10Fe10 | $306$ |  Mo40W20Pd10Gd10B10Cr10 | $5.53$
|  W30B20Au20Be10Nb10Ag10 | $296$ |  Ta20Nb20Al20Au20B10W10 | $5.06$
7 |  W20Co20Nb20Ag10B10Be10Mg10 | $284$ |  W20B20Ag20Nb10Si10Co10Pd10 | $5.24$
|  Cr20Ag20Ti20B10Gd10Be10Mg10 | $234$ |  Cr20Fe20W20Ca10B10Sn10Be10 | $3.78$

In Table 3, we list $10$ binary and ternary amorphous metal alloys selected
from $50\,000$ alloys considered in this study. By simple comparison one can
reveal that the predicted $E$ is mainly correlated with the mechanical
properties of chemical element with the highest mass fraction. For example,
for amorphous Cr80B20, we find $E\approx 305$ GPa, while the Young’s modulus
of pure crystalline Cr is $E\approx 279$ GPa. In the case of amorphous
W40Mo40B20, we have $E\approx 318$ GPa, while $E$ is $\sim 410$ GPa for pure
crystalline W. The mechanical properties can vary depending on the
concentration of the doped chemical elements and on the class to which these
elements belong (metals, nonmetals, lanthanides, etc.). For example, the
predicted Young’s modulus for Ni40Cr40Co20 is $E\approx 58$ GPa. At the same
time, the presence of Zr and Si in the Ni-based alloy doubles the Young’s
modulus $E$ (i.e. one has $E\approx 108$ GPa for Ni40Zr40Si20). For
Ni40Mo40W20, the machine learning model predicts $E\approx 183$ GPa, where the
refractory metals Mo and W are included [see Table 3]. The doping with
refractory metals, nonmetals and lanthanides (e.g. B, Si, Gd, La) makes it
possible to increase the strength of these alloys, which is actively used in
modern metallurgy to produce heat-resistant alloys [71, 72]. This simple
quantitative analysis confirms that the properly selected composition and
physical properties of the main chemical elements are most important in
determining the alloys that best match the required mechanical properties.

Table 3: Young’s modulus $E$ predicted by machine learning model for binary and ternary amorphous metal alloys. Here, the average accuracy is $\sim 87$ %. Alloy |  $E$, GPa | Alloy |  $E$, GPa
---|---|---|---
Cr80B20 | $305$ | Cu80Mg20 | $60$
W40Mo40B20 | $318$ | Cu60Mo40 | $154$
Ni40Cr40Co20 | $58$ | W40Ag40B20 | $234$
Ni40Zr40Si20 | $108$ | Cr40B40Gd20 | $217$
Ni40Mo40W20 | $183$ | Cr40Nb40La20 | $196$

##  5 Conclusions

In the present study, the machine learning model was applied to predict the
Young’s modulus $E$ and the yield strength $\sigma_{y}$ of amorphous metal
alloys with different compositions. More than $50\,000$ different alloys were
determined as well as $E$ and $\sigma_{y}$ were evaluated for each of them. It
was found that the artificial neural network trained on the basis of
information about the atomic number of a chemical element, its atomic mass,
covalent radius, ionization energy, electronegativity, thermal conductivity,
specific heat capacity, density, melting temperature and boiling temperature
allows us to correctly determine of $E$ and $\sigma_{y}$ of amorphous metal
alloys consisting $2$ to $7$ components and containing chemical elements with
atomic numbers from $Z=3$ to $Z=79$. Here, the mean relative error is
$\sim(12\pm 1)$% that is the good accuracy for the direct propagation
multilayer neural network. The results of the statistical treatment made it
possible to determine the chemical elements with the largest mass fraction,
whose presence in the alloy leads to a significant increase in the strength of
alloys. These chemical elements are B, Cr, Fe, Co, Ni, Nb, Mo, Pd and W. At
the same time, the quantities $E$ and $\sigma_{y}$ show a weak dependence on
the number of components in alloy. Thus, the most significant factors in the
synthesis of alloys with the desired mechanical properties are the properly
selected composition and the physical properties of the basic chemical element
of alloy.

## Acknowledgment

This research was funded by the Russian Science Foundation (project no.
19-12-00022).

## References

  * [1] Wang, W.H.; Dong, C.; Shek, C.H. Bulk metallic glasses. Mater. Sci. Eng. R Rep. 2004, 44, 45–89.
  * [2] Chen, M. A brief overview of bulk metallic glasses. NPG Asia Mater. 2011, 3, 82–90.
  * [3] Yeh, P.-Y.; Huang, J.C.; Jang, J.S.C.; Pan, C.-T.; Chen, C.-H.; Lin, C.-H. Recent Developments in Additive-Manufactured Intermetallic Compounds for Bio-Implant Applications. J. Med. Biol. Eng. 2022, 42, 800–815.
  * [4] Huang, L.J.; Lin, H.J.; Wang, H.; Ouyang, L.Z.; Zhu, M. Amorphous alloys for hydrogen storage. J. Alloys Compd. 2023, 941, 168945.
  * [5] Kruzic, J.J. Bulk Metallic Glasses as Structural Materials: A Review. Adv. Eng. Mater. 2016, 18, 1308–1331.
  * [6] Louzguine-Luzgin, D.V.; Polkin, V.I. Properties of bulk metallic glasses. Russ. J. Non-Ferr. Met. 2017, 58, 80–92.
  * [7] Li, H.X.; Lu, Z.C.; Wang, S.L.; Wu, Y.; Lu, Z.P. Fe-based bulk metallic glasses: Glass formation, fabrication, properties and applications. Prog. Mater. Sci. 2019, 103, 235–318.
  * [8] Malygin, G.A. Strength and plasticity of nanocrystalline materials and nanosized crystals. Phys. Usp. 2011, 54, 1091–1116.
  * [9] Galimzyanov, B.N.; Doronina, M.A.; Mokshin, A.V. Excellent glass former Ni62Nb38 crystallizing under combined shear and ultra-high pressure. J. Non-Cryst. Solids 2021, 572, 121102.
  * [10] Anikeev, S.G.; Artyukhova, N.V.; Shabalina, A.V.; Kulinich, S.A.; Hodorenko, V.N.; Kaftaranova, M.I.; Promakhov, V.V.; Gunter, V.E. Preparation of porous TiNi-Ti alloy by diffusion sintering method and study of its composition, structure and martensitic transformations. J. Alloys Compd. 2022, 900, 163559.
  * [11] Galimzyanov, B.N.; Mokshin, A.V. Mechanical response of mesoporous amorphous NiTi alloy to external deformations. Int. J. Solids Struct. 2021, 224, 111047.
  * [12] Lesz, S. Effect of cooling rates on the structure, density and micro-indentation behavior of the Fe, Co-based bulk metallic glass. Mater. Charact. 2017, 124, 97–106.
  * [13] Louzguine-Luzgin, D.V. Structural Changes in Metallic Glass-Forming Liquids on Cooling and Subsequent Vitrification in Relationship with Their Properties. Materials 2022, 15, 7285.
  * [14] Tournier, R.F.; Ojovan, M.I. Multiple Melting Temperatures in Glass-Forming Melts. Sustainability 2022, 14, 2351.
  * [15] Ojovan, M.I.; Tournier, R.F. On Structural Rearrangements Near the Glass Transition Temperature in Amorphous Silica. Materials 2021, 14, 5235.
  * [16] Galimzyanov, B.N.; Yarullin, D.T.; Mokshin, A.V. Kinetics of inherent processes counteracting crystallization in supercooled monatomic liquid. J. Phys. Condens. Matter 2022, 34, 454002.
  * [17] Kramer, J. Der amorphe Zustand der Metalle. Z. Phys. 1937, 106, 675–691.
  * [18] Brenner, A.; Couch, D.E.; Williams, E.K. Electrodeposition of alloys of phosphorus with nickel or cobalt. J. Res. Nat. Bur. Stand. 1950, 44, 109–122.
  * [19] Duwez, P. Structure and Properties of Alloys Rapidly. Quenched from the Liquid State. Trans. Am. Soc. Metals 1967, 60, 607–633.
  * [20] Pond Jr., R.; Maddin, R. A method of producting rapidly solidified filamentary castings. Trans. TMS-AIME 1969, 245, 2475.
  * [21] Peker, A.; Johnson, W.L. A highly processable metallic glass: Zr41.2Ti13.8Cu12.5Ni10.0Be22.5. Applied Physics Letters 1993, 63, 2342–2344.
  * [22] Demetriou, M.D.; Kaltenboeck, G.; Suh, J.-Y.; Garrett, G.; Floyd, M.; Crewdson, C.; Hofmann, D.C; Kozachkov, H.; Wiest, A.; Schramm, J.P.; Johnson, W.L. Glassy steel optimized for glass-forming ability and toughness. Applied Physics Letters 2009, 95, 041907.
  * [23] Conner, R.D.; Li, Y.; Nix, W.D.; Johnson, W.L. Shear band spacing under bending of Zr-based metallic glass plates. Acta Materialia 2004, 52, 2429–2434.
  * [24] Inoue, A.; Shen, B.L.; Chang, C.T. Fe- and Co-based bulk glassy alloys with ultrahigh strength of over 4000 MPa. Intermetallics 2006, 14, 936–944.
  * [25] Zhang, H.W.; Subhash, G.; Jing, X.N.; Kecskes, L.J.; Dowding, R.J. Evaluation of hardness–yield strength relationships for bulk metallic g lasses. Philosophical Magazine Letters 2006, 86, 333–345.
  * [26] Yuan, C.C.; Xi, X.K. On the correlation of Young’s modulus and the fracture strength of metallic glasses. Journal of Applied Physics 2011, 109, 033515.
  * [27] Mokshin, A.V.; Mokshin, V.V.; Mirziyarova, D.A. Formation of Regression Model for Analysis of Complex Systems Using Methodology of Genetic Algorithms. Nonlinear Phenom. Complex Syst. 2020, 23, 317–326.
  * [28] Mokshin, A.V.; Khabibullin, R.A. Is there a one-to-one correspondence between interparticle interactions and physical properties of liquid? Phys. A: Stat. Mech. Appl. 2022, 128297.
  * [29] White, A.A. Big data are shaping the future of materials science. MRS Bull. 2013, 38, 594–595.
  * [30] Rodrigues Jr, J.F.; Florea, L.; de Oliveira, M.C.F.; Diamond, D.; Oliveira Jr, O.N. Big data and machine learning for materials science. Discov. Mater. 2021, 1, 12.
  * [31] Agrawal, A.; Choudhary, A. Perspective: Materials informatics and big data: Realization of the “fourth paradigm” of science in materials science. APL Mater. 2016, 4, 053208.
  * [32] Tripathi, M.K.; Kumar, R.; Tripathi, R. Big-data driven approaches in materials science: A survey. Mater. Today: Proc. 2020, 26, 1245–1249.
  * [33] Shokrollahi, Y.; Dong, P.; Gamage, P.T.; Patrawalla, N.; Kishore, V.; Mozafari, H.; Gu, L. Finite Element-Based Machine Learning Model for Predicting the Mechanical Properties of Composite Hydrogels. Appl. Sci. 2022, 12, 10835.
  * [34] Merayo, D.; Rodríquez-Prieto, A.; Camacho, A.M. Prediction of Mechanical Properties by Artificial Neural Networks to Characterize the Plastic Behavior of Aluminum Alloys. Materials 2020, 13, 5227.
  * [35] Ciupan, E.; Ciupan, M.; Jucan, D.C. Determining the Mechanical Properties of a New Composite Material using Artificial Neural Networks. TJETT 2018, 66, 103–108.
  * [36] Galimzyanov, B.N.; Doronina, M.A.; Mokshin, A.V. Arrhenius Crossover Temperature of Glass-Forming Liquids Predicted by an Artificial Neural Network. Materials 2023, 16, 1127.
  * [37] Xiong, J.; Shi, S.-Q.; Zhang, T.-Y. A machine-learning approach to predicting and understanding the properties of amorphous metallic alloys. Materials & Design 2020, 187, 108378.
  * [38] Khakurel, H.; Taufique, M.F.N.; Roy, A.; Balasubramanian, G.; Ouyang, G.; Cui, J.; Johnson, D.D.; Devanathan, R. Machine learning assisted prediction of the Young’s modulus of compositionally complex alloys. Sci. Rep. 2021, 11, 17149.
  * [39] Amigo, N.; Palominos, S.; Valencia, F.J. Machine learning modeling for the prediction of plastic properties in metallic glasses. Sci. Rep. 2023, 13, 348.
  * [40] Xiong, J.; Zhang, T.Y.; Shi, S.Q. Machine learning prediction of elastic properties and glass-forming ability of bulk metallic glasses. MRS Commun. 2019, 9, 576–585.
  * [41] Galimzyanov, B.N.; Doronina, M.A.; Mokshin, A.V. Machine learning-based prediction of elastic properties of amorphous metal alloys. Phys. A: Stat. Mech. Appl. 2023, 617, 128678.
  * [42] Johnson, W.L.; Samwer, K. A Universal Criterion for Plastic Yielding of Metallic Glasses with a $(T/T_{g})^{2/3}$ Temperature Dependence. Phys. Rev. Lett. 2005, 95, 195501.
  * [43] Jeon, J.; Seo, N.; Kim, H.-J.; Lee, M.-H.; Lim, H.-K.; Son, S.B.; Lee, S.-J. Inverse Design of Fe-Based Bulk Metallic Glasses Using Machine Learning. Metals 2021, 11, 729.
  * [44] Inone, A. Bulk Amorphous Alloys: Preparation and Fundamental Characteristics (Materials Science Foundations). Trans Tech Publications, Netherlands 1998.
  * [45] Material properties database of the “Information Technologies in Physical Materials Science” Laboratory of Kazan Federal University. Available online: http://itphyms.int.kpfu.ru/ (accessed on 27.03.2023, by subscription).
  * [46] Scientific database of materials and their properties from Lawrence Berkeley National Laboratory. Available online: https://materialsproject.org/ (accessed on 27.03.2023).
  * [47] Wang, W.H. Correlations between elastic moduli and properties in bulk metallic glasses. J. Appl. Phys. 2006, 99, 093506.
  * [48] Wang, W.H. The elastic properties, elastic models and elastic perspectives of metallic glasses. Prog. Mater. Sci. 2012, 57, 487–656.
  * [49] Chen, X.-Q.; Niu, H.; Li, D.; Li, Y. Modeling hardness of polycrystalline materials and bulk metallic glasses. Intermetallics 2011, 19, 1275–1281.
  * [50] Qu, R.T.; Liu, Z.Q.; Wang, R.F.; Zhang, Z.F. Yield strength and yield strain of metallic glasses and their correlations with glass transition temperature. J. Alloys Compd. 2015, 637, 44–54.
  * [51] Chumachenko, K.; Iosifidis, A.; Gabbouj, M. Feedforward neural networks initialization based on discriminant learning. Neural Netw. 2022, 146, 220–229.
  * [52] Sharma, S.; Sharma, S.; Athaiya, A. Activation functions in neural networks. IJEAST 2020, 4, 310–316.
  * [53] Haykin, S. Neural networks and learning machines, 3rd ed.; Pearson Education, Inc., Upper Saddle River, New Jersey, USA, 2009.
  * [54] Li, J.; Cheng, J.; Shi, J.; Huang, F. Brief Introduction of Back Propagation (BP) Neural Network Algorithm and Its Improvement. Adv. Comput. Sci. Inf. Eng. 2012, 169, 553–558.
  * [55] Yan, Z.; Zhong, S.; Lin, L.; Cui, Z. Adaptive Levenberg-Marquardt Algorithm: A New Optimization Strategy for Levenberg-Marquardt Neural Networks. Mathematics 2021, 9, 2176.
  * [56] Sobhanifar, N.; Ahmadloo, E.; Azizi, S. Prediction of Two-Phase Heat Transfer Coefficients in a Horizontal Pipe for Different Inclined Positions With Artificial Neural Networks. J. Heat Transfer 2015, 137, 061009.
  * [57] $\breve{\text{S}}$ter, B. Selective Recurrent Neural Network. Neural Process. Lett. 2012, 38, 1–15.
  * [58] Turco, C.; Funari, M.F.; Teixeira, E.; Mateus, R. Artificial Neural Networks to Predict the Mechanical Properties of Natural Fibre-Reinforced Compressed Earth Blocks (CEBs). Fibers 2021, 9, 78.
  * [59] Marques, A.E.; Prates, P.A.; Pereira, A.F.G.; Oliveira, M.C.; Fernandes, J.V.; Ribeiro, B.M. Performance Comparison of Parametric and Non-Parametric Regression Models for Uncertainty Analysis of Sheet Metal Forming Processes. Metals 2020, 10, 457.
  * [60] Jetly, V.; Chaudhury, B. Extracting electron scattering cross sections from swarm data using deep neural networks. Mach. Learn.: Sci. Technol. 2021, 2, 035025.
  * [61] Malyshev, V.; Makasheva, A. The Relation between the Heat of Melting Point, Boiling Point, and the Activation Energy of Self-Diffusion in Accordance with the Concept of Randomized Particles. Open J. Phys. Chem. 2014, 4, 166–172.
  * [62] Kanematsu, N.; Inaniwa, T.; Nakao, M. Modeling of body tissues for Monte Carlo simulation of radiotherapy treatments planned with conventional x-ray CT systems. Phys. Med. Biol. 2016, 61, 5037–5050.
  * [63] Karpechev, E.V.; Pshenichnov, I.A.; Karavicheva, T.L.; Kurepin, A.B.; Golubeva, M.B.; Guber, F.F.; Maevskaya, A.I.; Reshetin, A.I.; Tiflov, V.V.; Topilskaya, N.S.; Cortese, P.; Dellacasa, G.; Arnaldi, R.; De Marco, N.; Ferretti, A.; Gallio, M.; Musso, A.; Oppedisano, C.; Piccotti, A.; Scomparin, E.; Vercellin, E.; Cical$\grave{\text{o}}$, C.; Puddu, G.; Siddi, E.; Szymanski, P.; Efthymiopoulos, I. Emission of forward neutrons by 158A GeV indium nuclei in collisions with Al, Cu, Sn and Pb. Nucl. Phys. A 2014, 921, 60–84.
  * [64] Agmon, N. Covalent radii from ionization energies of isoelectronic series. Chem. Phys. Lett. 2014, 595–596, 214–219.
  * [65] Matovi, B.; Yano, T. Silicon Carbide and Other Carbides. Handbook Adv. Ceram. 2013, 225–244.
  * [66] Braithwaite, E.R.; Haber, J. Molybdenum: An Outline of its Chemistry and Uses. Elsevier Science B. V., Amsterdam, Netherlands 1994.
  * [67] Barnhart, J. Occurrences, uses and properties of chromium. Regulatory Toxicology and Pharmacology 1997, 26, S3–S7.
  * [68] Rieth, M.; Dudarev, S.L.; Gonzalez de Vicente, S.M. et al. Recent progress in research on tungsten materials for nuclear applications in Europe. J. Nucl. Mater. 2013, 432, 482–500.
  * [69] Yamanaka, K.; Mori, M.; Kuramoto, K.; Chiba, A. Development of new Co–Cr–W-based biomedical alloys: Effects of microalloying and thermomechanical processing on microstructures and mechanical properties. Materials & Design 2014, 55, 987–998.
  * [70] Baldinozzi, G.; Pontikis, V. Phenomenological potentials for the refractory metals Cr, Mo and W. Journal of Physics: Condensed Matter 2022, 34, 315702.
  * [71] Sanin, V.V.; Kaplansky, Y.Yu.; Aheiev, M.I.; Levashov, E.A.; Petrzhik, M.I.; Bychkova, M.Ya.; Samokhin, A.V.; Fadeev, A.A.; Sanin, V.N. Structure and Properties of Heat-Resistant Alloys NiAl-Cr-Co-X (X = La, Mo, Zr, Ta, Re) and Fabrication of Powders for Additive Manufacturing. Materials 2021, 14, 3144.
  * [72] Wang, L.; Shen, J. Effect of heat treatment on the microstructure and elevated temperature tensile property of Fe-doped NiAl–Cr(Mo)–(Hf,Dy) eutectic alloy. Mater. Sci. Eng. A 2016, 654, 177–183.