Parameter estimation and fractional derivatives of dengue transmission model

1.   Introduction

In the early 1980s, Robert F. Furchgott and John V. Zawadzki, while analyzing the contractions and relaxations of vascular rings in response to certain pharmacological substances (adrenaline for constriction and acetylcholine for dilation), observed that acetylcholine did not produce vasodilation when the inner layer of cells in blood vessels was removed ^[1]. This relaxation or vasodilation, dependent on the vascular endothelium, was due to the release of a substance by acetylcholine that, when it reached the smooth muscles of the vascular wall, produced vasodilation. This unknown substance was called the endothelium-derived relaxing factor (EDRF); this substance had a very unstable character (in normal physiological conditions, 3 seconds of half-life), and the vasodilation that it produced was a consequence from activating a localized enzyme in the vascular wall called guanylate cyclase ^[2].

Identifying EDRF with nitric oxide (NO) was a major milestone in vascular physiology and pathophysiology research. Its discovery was the product of parallel work from different research groups that reached similar results ^[3,4,5,6], and it was at the end of the 1980s when researchers David S. Bredt, Solomon H. Snyder and Salvador Moncada began to relate NO to the nervous system in its role as a neurotransmitter ^[7,8,9].

NO is an unstable molecule that is free radical, gaseous and fat soluble and it is capable of crossing cell membranes without the aid of specific transporters or mechanisms. It is small ($1.15\ \ddot{A} = 1.15 \cdot 10^{-4} \mu m$), spreads rapidly and has a half-life between 3 and 5 seconds. It is produced from the amino acid L-arginine and its production is regulated by the enzyme nitric oxide synthase (NOS).

The two classes of NOS that exist, constitutive (cNOS) and inducible (iNOS), are distributed in very different and varied areas. Thus, the endothelial isoform (eNOS) of cNOS is located in the nerve fibres that surround blood vessels and endothelial cells, and the neuronal isoform (nNOS) is located in neurons of the central nervous system (CNS) and the peripheral nervous system (PNS). On the other hand, iNOS appears in glial cells and is only expressed in the event of brain injury or neurodegenerative diseases ^[10]. Based on the above, we observed that there are multiple neurons capable of synthesising NO, which is the subject of our modeling and simulation study. This NO acts in the CNS as an atypical neurotransmitter because it can be released by any part of the cell membrane without the need for presynaptic or postsynaptic structures, storage vesicles or transporter proteins.

The dynamics of NO are formed by three fundamental processes (generation, diffusion and self-regulation and recombination) ^[11] and begin through its generation, which is a process that occurs within the chemical communication itself and requires the existence of the NOS enzyme and $Ca^{2+}$ flux. This process causes a rapid and transitory release of moderate amounts of NO that diffuse in all directions. The diffusion of NO is governed by the gradient of its own concentration, and its objectives are those cells that, being within its reach, have the ability to recombine with it, such as the enzyme soluble guanylate cyclase (GCs) ^[12,13]. The maximum range of NO influence through diffusion can reach 300 $\mu m$, totaling approximately $2 \cdot 10^{6}$ synapses ^[14], and as this diffusion takes place, a process of self-regulation and recombination occurs with different substrates ^[15], with no reuptake processes. The generated NO contributes to inhibiting the activity of NOS by a negative feedback mechanism, leaving the area in a refractive period, during which NO generation does not occur again.

The presence of a molecule such as NO in the nervous system opens up new perspectives in the study of its functioning due to its various functions in both the PNS and CNS ^[11], identifying it as a possible underlying element and supporting the high capacity of adaptation, flexibility and prediction that the brain possesses. Confirming the latter is to understand the functioning of NO as a signaling molecule in the brain. The absence of decisive experimental data for this function has generated the development of various models. Thus, Tadeusz Malinski et al. ^[16] directly applied the Fick equation to model the behavior of NO through isotropic and homogeneous diffusion. Lancaster J. R. ^[17] associated the dynamics of NO with a random walk, intrinsically related to the parabolic equation of diffusion, and introduced the compartment concept to model NO. Wood J. and Garthwaite J. ^[18] based their model on an analytical solution of the diffusion equation and studied the physiological sphere of action and influence of an isolated source of NO, also analyzing the self-regulation process of NO. Vaughn M. W. et al. ^[19] characterized the dynamics of NO by means of a model composed of three regions with differentiated dynamics and followed the anatomical structure of the vascular endothelium (luminal, endothelial and adluminal). Philippides A. et al. ^[20] modeled the diffusion of NO from an irregular structure using an analytical solution of the diffusion equation and introduced a global morphology for a generation process where the generation of NO takes place in a spherical section. Suárez Araujo C. P. et al. ^[21,22,23], developed several models to study and determine the dynamics of NO from different approaches: An analytical model in a continuum based on the Fick equation that allowed working with cylindrical morphologies; a compartmental model that, using phenomenological equations of transport and compartmental analysis as an underlying formal mechanism, allowed considering different morphologies and characteristics of NO dynamics, such as the nonisotropy of the medium and the nonhomogeneous character of the basic processes of NO dynamics (generation and diffusion of NO and reaction of NO with other molecules).

Our model focuses on the dynamics of NO as a signaling molecule in VT, mainly on how NO transmits information using the extracellular environment in its three modes of operation: In the single synapse, synaptic spillover and Volume Transmission (VT) ^[24]. We do not rule out extending the scope of this model in the future, to incorporate the dynamics of NO recombination with other substances and other biochemical networks, including models that allow us to work from the perspective of the micro kinetics of NO.

In this field, metabolic models ^[25,26], represent a very important tool for our understanding of biological systems as they are more closely connected with cell behavior, and can adequately represent the effects of the dynamics of NO on target cells, which span different time scales and are frequently not recordable electrophysiologically ^[24].

Although there is a large body of work related to the dynamic of NO as a mediator of Central Nervous System (CNS) functions, the manner in which NO does this is not well understood, because its actions in target cells are mediated by metabotropic receptors with different time scales in its effects, and these are often not electrophysiologically detectable. On the other hand, we seek to know how a coherent and interpretable information transfer can be encoded by this signaling system in the brain. Particularly because the message itself is not channelled to any target, but is free to diffuse from its source in all three dimensions.

The main motivation of our work is to improve the knowledge and understanding the implications of NO in VT, and its effects in the formation of complex functional structures on which the VT performed by NO is supported, which are necessities for the synchronous functional recruitment of neuronal populations.

In this work, we model the effect of NO dynamics in its environment, as well as its implications in the formation of complex functional structures on which the VT performed by NO is supported, which is necessary for the synchronous functional recruitment of neuronal populations. Therefore, we propose a model based on the different changes that occur in the environment where this dynamic takes place. We use automata networks, which provide an alternative to classic and traditional models based on continuous dynamic systems-ordinary differential equations and partial derivative equations. The first model based on automata networks is deterministic, which presents a finite time of behavior based on an initial configuration of compartments that are in a state of NO generation.

This deterministic model is extendable and can incorporate complexity into its processes. One of these approaches allows the generation processes to occur arbitrarily in the compartments and be kept random over time. In modeling, fuzzy dynamics are considered in the behavior of the automata that are associated with these compartments. Thus, we obtain a model of a fuzzy automata network for the diffusion of NO.

We organize this work into two major sections, in addition to the introduction and the conclusions. Section 2 presents the formal tools and conceptualization of the automata networks (deterministic and fuzzy) used to model the NO dynamics. The results of the model, an analysis and study of deterministic convergence and an analysis of its dynamics are detailed in Section 3. We end this work with a section that summarizes the conclusions and identifies the methodological axis of our future work. The list of abbreviations and symbols is given in Table A1.1 (Annex 1. List of abbreviations and symbols).

2.   Method

The method of our work is based on modeling using the discrete mathematical structure of an automata network. In general, an automata network (AN) can be defined as a set of locally interconnected automata that evolve in discrete time steps through mutual interactions between them ^[27]. From a mathematical point of view, ANs are discrete dynamic systems.

Before going into the formal details of AN (in their two forms: deterministic automata networks and fuzzy automata networks), we define the basic components that compose them: Deterministic automata (DA) and fuzzy automata (FA).

Deterministic Automata (DA). A DA is a mathematical structure of states and transitions. The states of a DA represent a configuration in which a modeled real system is found with one or more states designated as starting state/s that represent the initial configuration of said real system.

Transitions in a DA are associated with changes in the configuration of the real system due to an action that can be internal or external. The former represent internal computing steps ($\tau$) and are not visible in the DA environment. An external action is visible in the DA environment and is used to interact with it.

Formally, a DA, denoted by A, is made up of the following four components:

$(ⅰ)$ A set of states $S_A$.

$(ⅱ)$ A nonempty set of starting states $I_A \subseteq S_A$.

$(ⅲ)$ Two sets of actions $V_A$ and $W_A$, external (and that interact with the DA environment) and internal, respectively, such that $V_A \cap W_A = \varnothing$ and where $Act_A = V_A \cup W_A$ defines the full set of actions associated with DA transitions.

$(ⅳ)$ A Cartesian product $\Delta_A \subseteq S_A \times Act_A \times S_A$, which defines the relationship of DA transitions.

Thus, we say that $s \overset{a}{\rightarrow } s'$, if $\exists (s, a, s') \in \Delta_A$, meaning that action $a$ is enabled in state $s$, and therefore there is a transition labelled $a$ from $s$ towards another state of DA.

The DA constitute a mathematical framework for the specification and analysis of real systems, but there are various formal extensions of them that collect aspects such as concurrency, asynchronous events, probabilistic configuration changes and environments with uncertainty. One of these extensions is probabilistic automata (PA), which has been formalized and developed by various authors ^{[27,28,29,30,31]}.

An additional extension to the above is fuzzy automata (FA), which are based on the concepts associated with fuzzy sets ^[32] and have been formalized and characterized by various authors ^{[33,34,35,36,37,38,39]}. DA and PA are particular cases of FA ^[34].

Fuzzy automata (FA). A FA is a DA where the final existence of the transitions between one state and the subsequent set of states depends on a fuzzy affiliation function that defines the degrees of transition. Therefore, in a FA, the transition from a state $s$ to a state $s'$, can be implemented if the value of the fuzzy function affiliation, before a given input $e$, meets a certain established existence criterion. Thus, a transition in a FA relates a state and an action to a fuzzy interval $[0, 1]$.

In this context, we define a fuzzy function affiliation (FFA) that characterizes a set $A$, which we also call a fuzzy set (FS), $\forall x \in X$, where we can calculate a value $f_A(x) \in [0, 1]$ and whose meaning is the degree of membership, or affiliation, that point $x$ has to the set $A$ and where $f_A(x) = 0$ corresponds to the fact that $x \notin A$, and $f_A(x) = 1$ with $x \in A$. Although the above can be interpreted as the simplest way to define membership for $\forall x \in X$ to FS $A$, there are other alternatives in which said membership is defined based on the parameters $\alpha$ and $\beta$, where $\alpha, \beta \in [0, 1]$, and $\alpha > \beta$. In this case, membership is defined based on the following three cases: 1) $x \in A$ if $f_A(x) \geq \alpha$, 2) $x \notin A$ if $f_A(x) \leq \beta$, and 3) the membership of point $x$ to FS $A$ is indeterminate when $\beta \leq f_A(x) \leq \alpha$.

Formally, we define an FA as an algebraic structure $A$, formed by the following five components:

$(ⅰ)$ A nonempty set of input states $E_A$.

$(ⅱ)$ A nonempty set of internal states $S_A$.

$(ⅲ)$ A nonempty set of output states $V_A$.

$(ⅳ)$ A first FFA, $f_A: S_A \times E_A \times S_A \rightarrow [0, 1]$ called a direct fuzzy transition function, where $f_A(s, e_j, s')$ defines the degree of state transition $s$ to the state $s'$, when the input $e_j$ is received.

$(ⅴ)$ A second FFA, $g_A: V_A \times E_A \times S_A \rightarrow [0, 1]$ called a direct fuzzy output function, where $g_A(v_i, e_j, s)$ defines the degree to which the automaton produces the output $v_i$ when it is in a state $s$ and receives an input $e_j$.

Based on the previous definition, we say that the transition $s \overset{a, f_A(s, e_j, s'), g_A(v_i, e_j, s)}{\rightarrow } s'$ occurs if $f_A(s, e_j, s')$ and $g_A(v_i, e_j, s')$ meet the membership criteria on the FSs $S_A \times E_A \times S_A$ and $V_A \times E_A \times S_A$.

The generalization of the FA behavior for an input sequence $E_k$ of arbitrary length $k$ can be defined by applying the definition of FFA composition ^[32].

In this case, $f_A(s, E_k, s') = f_A(s, e_0, s^1; s^1, e_1, s^{k-1}; …; s^{k-2}, e_{k-1}, s^{k-1})$ $= Max_{\{s^1, s^2, …, s^{k-2}\}} \cdot Min[f_A(s, e_0, s^1), f_A(s^1, e_1, s^2), …, f_A(s^{k-2}, e_{k-1}, s^{k-1})]$, obtaining in this case the degree of transition of the state s to the state $s^{k-1}$, when the sequence of inputs defined by $E_k$ is received. It is also necessary to apply the same process with the FFA $g_A$.

Automata Networks. Generalizing the above definitions associated with a DA and FA, an AN can be described as a function defined between two spaces, ${\bf{F}}: E \times S^{n} \rightarrow S^{n}$, where $E$ and $S$ are finite spaces called spaces of inputs and states. The AN is then formed on n interconnected automata, where the connection structure is defined by ${\bf{F}}$ in the following way: automaton $i$ receives a connection from $j$ if $F_i$ depends on variable $j$, where $F_i$ corresponds to component $i$ of ${\bf{F}}$.

The way in which automata are interconnected in the AN is given by their structure or topology of connections.

A state of the network is a vector ${\bf{x}}$ in $S^n$. The dynamics of the network are then defined as the rule that transforms said vector ${\bf{x}}$ of $S^n$ into a vector $y$ of $S^n$ when a given input ${\bf{e}}$ is received. The rule of parallel iteration is defined by ${\bf{y}} = {\bf{F}}({\bf{x}}, {\bf{e}})$ and can be interpreted as follows: In each time step, each automaton calculates its next state by means of the function $F_i$ on the current state ${\bf{x}}$ of the network. The AN iteration rule can present different modes of iteration: Parallel (the order in which each $F_i$ is calculated is irrelevant) and sequential (the order in which $F_i$ is calculated does matter).

As soon as the dynamics-or temporal evolutions-of the network are defined, the problem of its asymptotic behavior arises. Since the space of states S, on which the AN has been defined, is finite, all trajectories of the AN are periodic, i.e., from cyclical or fixed points.

A fuzzy automata network (FAN) is made up of n interconnected FA, where the connection structure is defined by the inputs that each FA receives.

Once the formal tools that we use in our modelling have been presented, we will define a series of concepts associated with NO dynamics and their effects on the diffusion environment.

Compartment. The concept of compartment determines, through a mathematical construct, the minimum place study where there is a complete NO dynamic diffusion defined by the automaton ^[22,23], which is part of the neuronal substrate. From a biochemical point of view, it could be identified as a chemical computing environment capable of producing the different processes involved in NO dynamics: NO generation, NO reception and NO self-regulation ^[11], as shown in Figure 1(a). In direct correspondence with each of these processes and to identify the set of possible states in which said compartment may exist, it is established that such processes may be in a state of activation or operation or in a state of nonactivation.

We can then match the situation in which the NO generation process is active in the compartment with that in which certain chemical machinery is active in the biological substrate that converges in the generation of NO. Likewise, we can say that when NO reception is active in a compartment, a combination of NO and its possible receptor molecules is occurring in the biological substrate associated with the compartment, and thus some functional and/or metabolic change is being activated with said NO reception, as shown in Figure 1(b).

We treat the self-regulation process independently, knowing that, at the biological level, NO is depleted as it diffuses and combines with different substrates and that certain levels of NO condition its own generation. We can see in the definition of the set of states and transition function that the state of self-regulation is always later than that of generation and that when a compartment is in that state, there is no NO. We can also identify how, at the level of the neuronal substrate, we will have a state of transmission, which is a direct consequence of the state of diffusion of NO that occurs in its dynamics. In this state, the compartment receives NO from the environment and generates NO that will be received in other compartments, Figure 1(b).

Neighborhood and Scope. Within our set of compartments, whose structural and functional definition is established in the following two sections, we can define the following relationship:

Relationship belonging to a neighborhood. Given two compartments $C_i$ and $C_j$, we say that $C_i$ belongs to the neighborhood of $C_j$ if and only if $C_i$ is influenced by the state of $C_j$. In this case, we say that $C_j$ creates the neighborhood $\Pi_{C_j}$ and that $C_i \in \Pi_{C_j}$.

The above relationship allows us to define the neighborhoods of each compartment, and these are the ones that inherently define the scope of NO that is generated in each compartment to later diffuse. In Figure 2, we can see how compartment $C_i$ simultaneously belongs to the neighborhoods $\Pi_{C_j}$, $\Pi_{C_k}$ and $\Pi_{C_p}$. This means that the dynamics of compartment $Ci$ are influenced by the states in which compartments $C_j$, $C_k$ and $C_p$ are found. Likewise, compartment $C_i$ will create a neighborhood $\Pi_{C_i}$ to which compartments $C_j$, $C_k$ and $C_p$ may belong.

2.1.   Deterministic automata network (DAN) for modeling the effect of NO in its diffusion environment - states and transition function

Based on the previous definitions, we identify the following states in which the automaton can be found:

(1) Nonactivity state, (n). A compartment is in a nonactive state when the concentration of NO, or its variation, is negligible.

(2) Receiving state, (r). A compartment in a receiving state is receiving NO. This NO is generated or transmitted by any of the compartments that are creating the neighborhoods to which it belongs.

(3) Transmission state, (t). Two different situations can cause a compartment to be in a transmitting state.

– The first of these situations occurs when the generation of NO occurs in a compartment. When (as we will see later) and if there is a compartment in its neighborhood in a receiving state, this would imply that our automaton would go to a state of NO generation. However, if it also happens that the generation of NO occurs in any of the compartments, to whose neighborhood this compartment in question belongs, then the state associated with this situation is that of transmission.

– The second of the situations presented corresponds to the opposite sequence. A compartment may be receiving NO and begin a process of NO generation in it, which takes it to the transmission state if there is a compartment in its neighborhood with the capacity to receive NO.

The state of transmission, in the definition of our automaton, has been called this way because it intrinsically captures the existence of a relationship between what is received and what is generated from NO. It is not known if the reception of NO can have some implication in subsequent generation processes. In the same region of the nervous tissue, we understand that this relationship may occur in higher level processes. Likewise, it is not known what NO reception in a certain region of nervous tissue may imply when there is already a process of NO generation in that region. Our hypothesis could be an acceleration towards self-regulation. Therefore, this name identifies only one state of our automaton that is reached by the two situations previously explained, and from the point of view of transmission of the information, independent collection is relevant.

Figure 3 shows a sequence of states with which to better understand the transmission state. A small section of an automaton network is shown where compartments $C_j$ and $C_p$ have been highlighted along with the neighborhoods $\Pi_{C_j}$ and $\Pi_{C_p}$ that they form and thus the scope that the NO generated in compartments $Cj$ and $Cp$ can have. In Figure 3(a), we can see the initial state of the compartments and their associated neighborhoods. Figure 3(b) shows how the activation of a NO generation process in compartment $C_j$ initiates a receiving state in all the compartments that belong to its neighborhood, including $C_p$ itself. Figure 3(c) shows how a generation process is activated in compartment $C_p$, thus initiating a transmission state (in the same way, compartment $C_j$ also goes to a transmission state), without any relationship existing between the generation that has occurred and the NO that it received previously. Figure 3(d) shows how the generation of NO in $C_j$ ends, and thus $C_j$ goes to a state of self-regulation and $C_p$ to one of generation.

(4) NO generation state, (g). A compartment is in this state when NO generation takes place in it. Through this process of NO generation, the compartment influences those compartments that belong to its neighborhood.

(5) State of self-regulation, (a). The state of self-regulation is subsequent in sequence to the state of NO generation. During this state, there are no dynamics related to NO. It is important to bear in mind that, both in the self-regulation sate and in the receiving state, NO disappears. The self-regulation of NO is a process by which the generated NO disappears. At this time, it is being argued that this disappearance is carried out in proportion to the amount of NO that exists in the environment, although some models assume this occurs with other types of dependence ^[19]. The automaton model that we present here does not handle amounts of NO, and self-regulation is introduced as a state after the generation and forced transition.

We begin the specification of the transition function of the automaton by defining the jump behavior between each of the possible states in which said automaton can be found.

Two versions of the transition function are proposed, as shown in Table 1 (version 0 and version 1). Both versions reach a finite and deterministic behavior regardless of the initial configuration, and the generation and nonactivity states of NO are possible only in these initial configurations. These versions, although they present a complex behavior, should be seen as an intermediate step for the final definition of our fuzzy transition function, which will require stochastic behaviors to reflect NO generation and self-regulation more realistically.

Figure 4 shows the state diagrams of both versions of the transition function, which cause the environment to be divided into three different basic dynamics: The dynamics of nonactivity or receiving NO, the dynamics of NO generation and the dynamics of NO transmission. The possibility of passing between the previous dynamics is the main difference between the two versions that have been defined as the transition function.

As seen in the definition of the transition functions in Table 1, in version 0 of our automata network (rules $r_0$, $r_1$ and $r_2$), whether a compartment is in the dynamic of nonactivity or NO reception, or in the dynamics of generation or transmission of NO, will depend mainly on the initial configuration, since in no rule of its transition function is it allowed that a compartment initiating the receiving dynamics can pass to the dynamics of generation or transmission. We can also see that the reverse transition is allowed in this version of the automaton. A compartment can start in a generation dynamic and later move to a receiving dynamic directly or by a temporary transition through the transmission dynamic.

Version 1 of our automata network (rules $s_0$, $s_1$, $s_2$ and $s_3$) completes the previous scenario, allowing the transition between the dynamics of nonactivity or reception, originating in an initial configuration, to generation and transmission dynamics.

Therefore, our two versions of the transition function have 7 rules (defined algebraically in Table 1) to control the previous behavior. Below is a detailed textual definition of these rules:

Version 0:

● Rule $r_0$: Represents a situation in which, being a compartment without NO dynamics, it begins to receive NO. The automaton associated with said compartment will change from a state of nonactivity (n) to a state of reception (r). This occurs when any of the generating automatons of the neighborhoods to which the associated automaton belongs are in a generation (g) or transmission (t) state. The automaton will continue in this state of reception (r) as long as the previous situation is maintained. When this situation ceases to occur, the automaton returns to a state of nonactivity (n).

The rule $r_0$ controls the dynamics of those zones in which NO is not generated and only NO is received. The areas where it is produced and received will NOT be controlled by Rules $r_1$ and $r_2$.

● Rule $r_1$: Rule $r_1$ controls passing from NO generation dynamics to transmission dynamics, where there is NO generation and reception.

Therefore, we are faced with a situation where some of the generating automatons of the neighborhoods to which the associated automaton belongs are in a generation (g) or transmission (t) state, and some of the automatons of their own neighborhood are found in a state of reception (r) or nonactivity (n).

When the above occurs, the automaton goes from a generation state (g), which it reached as a result of its initial configuration, to a transmission state (t), remaining in said transmission state (t) as long as said rule is met.

● Rule $r_2$: This rule identifies a situation after an NO generation dynamic, which in our model is associated with an NO self-regulation state. There is a period of nonactivity in the compartment, which is necessary after the dynamics of NO generation and transmission.

With the automaton in a generation state (g), we go to a self-regulation state (a) when the automaton is being influenced by states of generation (g) or transmission (t) from the automatons that have it in their neighborhoods, and in the neighborhood of the automaton itself, there are also automatons in a state of generation (g) or transmission (t).

When the conditions of Rules $r_1$ and $r_2$ are not met, the automaton must continue in a state of generation (g).

Version 1:

● Rule $s_0$: Represents a situation in which, being a compartment without NO dynamics, it begins to receive NO. Therefore, the automaton will change from a state of nonactivity (n) to a state of reception (r) and remain in this state as long as it continues to receive NO. This occurs when any of the generating automatons of the neighborhoods to which the associated automaton belongs are in the generation (g) or transmission (t) state, and none of the automatons that belong to its neighborhood are already in a reception state (r). When this situation ceases to occur, the automaton returns to a state of nonactivity (n).

● Rule $s_1$: This rule identifies the conditions that must be produced to cause a change between the dynamics of nonactivity or NO reception and the dynamics of NO generation or transmission.

With the automaton in a state of nonactivity (n), we go to a state of generation (g) when the automaton is being influenced by states of generation (g) or transmission (t) from the automatons that have it in its neighborhood, and there are automata from its own neighborhood in a state of reception (r), a situation in which there is a functional demand for NO.

● Rule $s_2$: Like Rule $r_1$ of version 0 in the transition function, this rule controls a transition from the dynamics of only NO generation to the dynamics where there is generation and reception (transmission).

We are faced with a situation where if any of the generating automatons of the neighborhoods to which the automaton belongs are in a generation state (g) and some of the automatons from its own neighborhood are in a reception state (r) or nonactivity state (n), the automaton goes to the transmission state (t).

If the above condition is not met, the automaton goes to a state of self-regulation (a).

● Rule $s_3$: This rule is directly associated with a condition where the compartment should leave the NO transmission state (t).

If any of the generating automatons of the neighborhoods to which the automaton belongs are in a state of transmission (t), and any of the automatons from its own neighborhood are in a state of reception (r) or nonactivity (n), the automaton goes from a transmitting state (t) to a receiving state (r).

If the above condition is not met, the automaton goes to a state of self-regulation (a).

3.   Results and discussion

Of the wide spectrum of particularities presented by NO dynamics, such as 1) nonlocalized generation, 2) existence of a maximum range of NO beyond which there should be no influence thereof, 3) existence of areas where NO is recombined with other substances to influence the dynamics of other mechanisms of cellular communication, and therefore, in higher order processes such as neural plasticity and learning and 4) existence of substance self-regulation, among others, in our work, we will focus on how NO performs complex structures formation when it develops its dynamics and on which the possible transmission of information is constituted, which is necessary to provoke the synchronous functional recruitment of involved neural populations.

It is in this section that the result of analyzing the modeling through the DANs made explicit in the previous sections can be observed. This will allow us to know what type of behavior they develop throughout the generations.

We focus on the 1D versions of our automaton networks for modeling NO, and these studies are totally valid for automata networks of different dimensionalities ^[40]. Likewise, these results will be extrapolated to their 2D versions.

Figure 5 shows different 1D evolutions associated with the two versions of our automaton network; Figure 5(a), (b) for version 0; Figure 5(c), (d) for version 1. In this case, networks with 32 automata are used. The initial configurations have been established randomly, following the premise that each automaton can be either in a generation state (g) or nonactivity state (n) with equal probability, assuming the above that in these initial configurations, we will have approximately 50% of the automatons in a generation state (g). In these figures, it is observed how, in all cases, separate structures of a stable or periodic type are generated.

In version 0 of the DAN, the structures are supported by the different convergence cycles of each automaton, and we can quantify their appearance based on the level of NO generation contained in the initial configurations, as well as on the range of the neighborhood. In Figure 6, we can observe the number of times that the convergence cycle appears, composed only of the generation state (g) and for different reaches of the neighborhood: Range 1 for a neighborhood with three cells of the lattice (the two neighbors and the cell itself), range 2 for 5 cells and range 3 for 7 cells. This figure shows us how the occurrence of this convergence cycle varies as a function of the probability of the generation in the initial configuration and the scope of the neighborhood. It is observed that an interval occurs in which said convergence cycle disappears completely, the latter being able to be interpreted as a control mechanism of the NO level.

Annex 2 (identification and analysis of the different convergence cycles of the automata networks for modeling NO) compiles this quantification for all the possible sequences of these convergence cycles of said version 0 for our automaton network.

As we have already indicated, in version 0 of the transition function for the automaton network, the generation state (g) of a compartment can only occur in the initial configuration of the automaton. This is the main difference between the two versions of the developed automaton network, since in version 1, the process of generating NO in the automaton cannot only occur in the initial configuration.

As seen in the regular expressions that define the sequences of states that can be produced,

Both expressions are mutually self-contained, and therefore, the automaton can begin in a state of nonactivity (n) and later go to a state of generation (g) depending on whether the rules that define its transition function are met. This makes a detailed analysis of all the possible sequences, or convergence cycles, not feasible, as presented by version 1 of our automaton network.

3.1.   Qualitative and quantitative analysis of automata

Once this first visual analysis of both versions of our automaton network has been carried out, a qualitative analysis is carried out, taking as a cataloguing reference the classification established by Stephen Wolfram ^[40]. In this classification, there are four different types of dynamics for automata networks, which, depending on their variation in space and time from a random initial configuration, are classified as follows:

● Class Ⅰ. The evolution of the automaton network converges to a homogeneous state, without spatial or temporal structures of any kind.

● Class Ⅱ. The evolution of the automaton network tends to separate structures of a stable or periodic type.

● Class Ⅲ. The evolution of the automaton network presents chaotic patterns. Fractal structures emerge spatially, and cycles of very long length are observed.

● Class Ⅳ. The evolution of the automaton network generates localized complex structures, which spread and whose duration increases exponentially with the size of the network.

The first three classes correspond qualitatively to the three types of behaviors observed in continuous systems (attractors, periodic/quasiperiodic and chaotic).

From the qualitative analysis, the temporal evolutions of both versions of our automaton network are determined to be Class Ⅱ.

Quantitative analysis is driven by the work of Chris G. Langton ^[41]. First, the set $D_N^K$ of all the possible transition functions follow the algebraic structure $\Delta : \Sigma^N \rightarrow \Sigma$, where $K$ corresponds to the number of states of the automaton and $N$ to the number of neighbors that are involved in said transition function, including the automaton for which the transition function calculates the new state.

Let $s_q$ be a state that we identify as a quiescent state. Once this quiescent state is identified, we can identify the $n$ transitions that lead us to that state. Once identified, we are interested in quantifying the values ($K^N - n$) because they represent the rest of the transitions. The degree of heterogeneity in the behavior of the automaton will depend on the way in which said ($K^N - n$) transitions are defined. For this purpose, the parameter $\lambda$ is established according to expression 3.2 ^[41]:

If $n = K^N$, all transitions have sq as the final state, and the behavior of the automaton is completely homogeneous (all initial configurations end in $s_q$) and $\lambda = 0$. On the other hand, if $n = 0$, no transition has $s_q$ as the final state, and then $\lambda = 1$. The most heterogeneous behavior occurs when all state transitions (including those with $s_q$ as the final state) are equally represented, which happens when $n = K^N-1$ and $\lambda = 1-1/K$. Based on the definition of the parameter $\lambda$, we perform its calculation for the two versions of the automaton network when we work in a 1D environment.

As seen in Tables 2 and 3, for each possible state $s_i(t)$ and states of its neighbors $s_{i-1}(t)$ and $s_{i+1}(t)$, the state $s_i(t+1)$ is made to correspond according to the rules that define our transitions, as shown in Table 1. In this case, we are working with a value of $K = 5$, corresponding directly with the states: Nonactivity (n), generation (g), transmission (t), reception (r) and self-regulation (a). The definition of the neighborhood of the automaton and the fact that we are working in a 1D environment give us a value of $N = 3$.

Based on the logical rules defined in Table 1, all possible transitions are detailed in Tables 2 and 3 as a preliminary step for calculating the parameter $\lambda$. The expression of the transition functions, according to these tables, is defined as follows: If $[s_(i-1)(t) s_(i+1)(t) = 'nr' ] \wedge [s_i(t) = 't']$ then $s_i(t+1) = 'a'$, following, in this case, the transition indicated in Table 3.

The above allows us to calculate the value of $\lambda$ for a certain quiescent state $s_q$ by direct application of the expression ($K^N$ - number of occurrences of $s_q$)/$K^N$, obtaining the values shown in Table 4 for both versions of our automata network.

Taking as $s_q$ the state of nonactivity, we have $\lambda_{v_0} = 0.656$ and a $\lambda_{v_1} = 0.632$, as shown in Table 4. These values allow us to locate the transition functions of our automata in the set of all possible functions, $D_N^K$. It is important to note that when any of the other states is selected as state $s_q$, the values of $\lambda_{v_0}$ and $\lambda_{v_1}$ change.

To analyze the emergence of complexity produced by the dynamics of our automata networks, we will use the probabilistic approach offered by entropy as a basic measure of self-information. For a discrete process $T$ of $K$ possible states, this is defined according to Eq (3.3).

where $p_i$ corresponds to the probability of state $i$ occurring in process $T$.

On the other hand, to quantify the degree of cooperation that may exist in our DANs, on which the synchronous functional recruitment feature is supported, we quantify the level of correlation that exists in the events that occur (state changes) in the automata. For this, the concept of mutual information $I(T_n, T_m)$ is used between two automata $n$ and $m$, in which the discrete processes $T_n$ and $T_m$ occur. The magnitude is defined as a function of the individual entropy $H(T_n)$ and $H(T_m)$ of the two automata and the entropy of the two automata considered as a joint process $H(T_{n, m})$.

Therefore, the mutual information is given by the following expression:

This measure will have direct dependence on the correlation of process $T_n$ with the state of process $T_m$. Thus, high values in the average of $I(T_n, T_m)$ will imply a high cooperation between automata $n$ and $m$. In contrast, a functional independence, or change of states, between automata will assume low values of the previous measure.

By virtue of the defined magnitudes, our network of automatons for NO modeling adequately incorporates the characteristics of complex structure formation by having an intermediate value of the average of the general entropy $\overline{H}$, understanding that the complex is between the order of a system, where $\overline{H} \approx 0$, and the total disorder, where $\overline{H}$ presents its highest values. On the other hand, having a high synchronous functional recruitment between the various dynamics of NO implies that we must have high values in the average general mutual information $\overline{I}$.

It should be taken into consideration that since the previous magnitudes are dependent on the trajectories followed by the automata network, according to the initial configuration, and since the latter has a random character, in our quantitative study, average values of these magnitudes are calculated throughout a set of executions where the initial configuration changes in its form for each one of them, but the proportions of automata that can be in an initial state of generation (g) or nonactivity (n).

In the same way, all the calculated variables that depend on the trajectories that the network of automatons follow are organized on the axis of $\lambda$ associated with the set $D_N^K$, being able to locate our networks of automatons in said set for comparison.

Figure 7 shows us the values of the general entropy $\overline{H}$ and general mutual information $\overline{I}$ for a subset of $D_N^K$ along the dimension $\lambda$ in the range between 0 and $1-1/K$. This interval corresponds to the interval that goes from the values of $\lambda$ associated with the most homogeneous automata networks $\lambda = 0$ to the most heterogeneous $\lambda = 1-1/K$.

Figure 7(a) shows how the value of $\overline{H}$ in version 1 of our automaton is approximately $\overline{H} = 1.25$, which is approximately half the value ($\overline{H} \approx 2.1$) that the automata networks of $D_N^K$ have for that same $\lambda$. This fact indicates that version 1 of our automata network forms complex structures in its dynamics.

Figure 7(b) shows the suitability of version 1 of our automata network from the perspective of synchronous functional recruitment, since this network has a value of $\overline{I} = 0.8$, which practically doubles the average value of the rest of the automata networks of $D_N^K$.

In both comparisons shown in Figure 7, for the selection of the rest of the automata networks that make up the subset of $D_N^K$, the "table-walk-through" procedure has been carried out ^[41], taking as a seed version 1 of the automaton network and modifying its transition function stochastically to obtain automaton networks with different $\lambda$ within the study interval indicated above, and a state $s_q$ equal to that of nonactivity (n).

In Figure 8, we have the same comparison of the values of $\overline{H}$ and $\overline{I}$, but in this case, the subset of $D_N^K$ has been generated with the "table-walk-through" process, using version 1 as a seed automaton network and the state $s_q$ = self-regulation (a). The main difference is if we compare the levels of $\overline{H}$ (Figures 7(a) and 8(a)), when we work with $s_q$ = self-regulation (a), its lowest values are close to that of our automata network, implying that the level of complex structures of our automata network does not present a differentiating character for the rest.

Figure 9 shows a relationship factor that has $\overline{H}$ and $\overline{I}$ with the temporal evolution of the automaton networks of $D_N^K$ for different values of $\lambda$. This figure shows the evolution of 6 automata networks made up of 32 automata in a cyclical lattice, where their transition functions have been obtained by taking both versions of our automata network for modeling NO and where $s_q$ corresponds to the state of nonactivity (n). Figure 9(a)–(c) correspond to version 0, and Figure 9(d)–(f) correspond to version 1. It can be seen in these figures how the temporal evolution of the states of these networks presents structures of a certain complexity and fractal tendency, and they are characterized by different values, not only for the parameter $\lambda$ but also with regard to the values of $\overline{H}$ and $\overline{I}$.

The behavior observed in the analyzed DANs is also present in larger automaton networks, as we can see in Annex 3. Detail of the evolution of the automata networks with 128 automata (version 0 and version 1) of the automata network for modeling NO, for DANs with 128 automata.

In the developed quantitative analysis, it has also been verified how dependent the values of $\overline{H}$ and $\overline{I}$ are in relation to the initial configurations and to the level of NO generation that may exist in them. Figure 10 shows the evolution of $\overline{H}$ and $\overline{I}$ when we vary the percentage of automata that are in a generation state (g) in the initial configuration.

In this figure, it can be seen that the average values of $\overline{H}$ show a slight increase (from a value of $\overline{H} \approx 1.2$, up to a value of $\overline{H} \approx 1.3$) as we increase the level of NO generation present in the initial configurations, from 10% to 70%, becoming unstable for percentages greater than the latter. For the case of $\overline{I}$, we see that its variation, despite being upwards, is practically negligible, staying around $\overline{I} \approx 0.8$ and behaving in the same way as $\overline{H}$ once the NO generation percentage exceeds 70%.

On the other hand, it seems that the storage of information by any system (be it discrete or continuous) implies a low entropy, and the transmission of information implies the existence of an increase ^[42]. Likewise, a high mutual information supposes a high correlation between the automata.

The final analysis for version 1 of our automaton network seeks to see what position the latter has in the plane facing both magnitudes: Entropy versus mutual information, which is always in compared with the rest of the automaton networks of $D_N^K$.

Figure 11 shows that version 1 of our network is located in the area of the plane that gives it a medium entropy level. Consequently, we have a high formation of complex structures, and on the other hand, the level of mutual information seems to be at high values compared to all the automata networks that are defined by $D_N^K$. The above allows us to argue that version 1 of the automaton network also presents high levels of synchronous functional recruitment.

A more exhaustive analysis of how version 1 of our automaton network behaves in relation to the values of $\overline{H}$ and $\overline{I}$ for all the possible initial configurations is shown in Figure 12. The different colors determine the level of NO generation that the initial configurations present. It is observed in this figure that the initial configurations that make our network of automatons exit the area given by the following values are minimal ($\overline{H} \approx 0.3$ and $\overline{I} \approx 0.8$), which is identified as suitable for having an adequate level of complex structure formation and synchronous functional recruitment. It is also seen in this figure that the initial configuration that produces maximum values of $\overline{H}$ and $\overline{I}$ corresponds to an initial configuration composed of the following sequence of states: nnggnngg… nngg, where n and g correspond to the nonactivity and generation states respectively and whose evolution is shown in Figure 13.

We extend our study to 2D DAN, where results parallel to those achieved with 1D DAN can be observed. Thus, Figure 14 shows the complex structure formation by all the dynamics associated with the different states, Figure 14(a)–(e), in which automata may exist when we work with a network of 10,000 automata arranged in a $100 \times 100$ lattice and where each automaton has 5 neighbors, in correspondence with a Moore-type neighborhood. In these figures, the formation of a set of structures is perceived when the automaton network is in the 75th generation.

Figure 15 shows the convergence of the same automata network with version 1 of the transition function towards the final complex structures when we sufficiently advance in the number of generations. The edges of the structure show a changing and cyclical behavior.

The implementation of the different versions of the automaton network was developed in Python 3.8.13 (www.python.org) using the TensorFlow 2.3.0 libraries (www.tensorflow.org). The source code used in all the experiments in this article is available for download on GitHub (https://github.com/pablo-fernandez-lopez/NANetwork_NO).

4.   Conclusions

In this work, a model based on deterministic automata networks is proposed for modeling the effect that NO dynamics exert on the environment through which it diffuses in its role as a molecule with the ability to perform VT. We carry out the formalization of its transition function through the logical extrapolation of those mechanisms associated with the diffusion dynamics of NO as a neuroactive substance. The obtained model adequately depicts the characteristics of complex structure formation and synchronous functional recruitment.

To achieve this, we have built two versions of transition functions (version 0 and version 1) that segment the environment into three different basic dynamics: the dynamics of nonactivity or NO reception, the dynamics of NO generation and the dynamics of NO transmission. The possibility of passing between the dynamics of nonactivity or NO reception to the rest of dynamics is the main difference in version 1 compared to version 0. The two versions of the transition function, defined and analysed in this work utilize a deterministic NO generation process. This generation of NO occurs in the initial configuration of the automata (case of version 0) or in the initial configuration by meeting a specific logical rule of its transition function that responds to the situation of the functional requirement of NO (case version 1).

Both versions cause separate structures of a stable or periodic type in all the sequences of states through which the automata network passes for each initial configuration. These structures are supported by the different convergence cycles that each automaton develops, whose occurrences depend on the level of NO generation of the initial configuration as well as the range of the neighborhood. From a qualitative perspective, the two versions of our automaton network are classified as class Ⅱ on the Stephen Wolfram rating scale ^[40].

The quantitative analysis of version 1 of our automaton network, when compared with the rest of the possible automata networks generated and organized according to the heterogeneity measurement parameter $\lambda$, defined by Chris G. Langton ^[41], presents adequate values for entropy and mutual information ($\overline{H} \approx 1.3$ and $\overline{I} \approx 0.8$), achieving an adequate predisposition of the network for complex structure formation and synchronous functional recruitment necessary to model the VT and study its implications in mechanisms and higher processes of the brain, such as learning and memory formation.

Working with version 1 of our automaton network in 2D environments, it is observed that NO dynamics produce areas of isolation and segmentation of the environment in relation to the characteristics of complex structure formation and synchronous functional recruitment. These complex structures present a zonal convergence when we sufficiently advance the number of generations, where the edges of the structure present a changing and cyclical behavior.

Finally, we propose the first discrete model in all its variables that can work with different NO dynamics and analyze the implications of VT in more complex architectures and aspects related to learning and memory formation.

We consider the DAN model proposal presented in this work, in its two versions, as a first step, and we identified the need to extend our model to incorporate stochastic conditions that make the state of NO generation be induced by higher mechanisms or brain processes. To carry out this generalization, which will also constitute a model that can accommodate arbitrary processes in decision-making mechanisms, we propose the use of fuzzy automata networks. This model will be part of a complete formal framework of volumetric transmission in the brain and in artificial neural networks and therefore in complex decision-making systems.

Use of AI tools declaration

The authors declare that they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The research presented in this paper has been funded by the Project "Investigación en Computación Neuronal por el Grupo de Investigación CIPERBIG (Research in Neural Computation by the CIPERBIG Research Group) (ULPGC)"; No: 23/2021, from Cabildo de Gran Canaria.

We are thankful to the "Council of First Vice-presidency and Public Works, Infrastructures, Transport and Mobility of the Cabildo de Gran Canaria".

Conflict of interest

The authors declare no conflicts of interest.

Supplementary

Annex 1.   List of abbreviations and symbols

Annex 2.   Identification and analysis of the different convergence cycles of the automata networks for modeling NO

Version 0 of the automata network:

Annex 3.   Detail of the evolution of the automata networks with 128 automata (version 0 and version 1) of the automata network for modeling NO