*2.3. The Model of a Network*

An example illustrating the idea of a considered computing oscillator network is shown in Figure 3. The network is formed by three coupled chemical oscillators (nodes) graphically represented by circles. The nodes have different functions in networks. The upper one (#1) is a normal node whose illumination time is fixed. Two bottom nodes (#2) and (#3) are inputs of *x* and *y* coordinates, respectively. The arrows interlinking the oscillators represent reactions that exchange the activators between nodes. The arrows directed away mark the activator decay.

**Figure 3.** The idea of a computing oscillator network. Circles represent network nodes that are chemical oscillators. The nodes have different characteristics. The upper one (#1) is a normal one, and its illumination function is fixed. The bottom nodes (#2) and (#3) are inputs of *x* and *y* coordinates, respectively. The arrows interlinking oscillators represent reactions that exchange the activators between nodes. The arrows directed away mark activator decay (reaction 3).

If the *j*th oscillator is a normal one, then the value of *tillum*(*j*) that appears in the definition of *φ*(*t*) (Equation (4)) does not depend on input values. Normal nodes are supposed to moderate the interactions within the network. The set of their illuminations can be regarded as a program executed by the network. If an oscillator is considered the input of a predictor *p<sup>i</sup>* . then the value of *tillum*(*j*) is functionally related to *p<sup>i</sup>* . For the analysis presented below, it is assumed that this function is defined by two parameters, *tstart* and *tend*, and has the form:

$$t\_{illum}(j) = t\_{start} + (t\_{end} - t\_{start}) \cdot p^i \tag{5}$$

Keeping in mind the symmetry of the Japanese flag problem, the values of *tstart* and *tend* are the same for all predictors corresponding to *x*- and *y*-coordinates. This means that for record (*xn*, *yn*,*sn*):

$$t\_{illum}(2) = t\_{start} + (t\_{end} - t\_{start}) \cdot x\_n \tag{6}$$

and

$$t\_{illum}(\mathfrak{J}) = t\_{start} + (t\_{end} - t\_{start}) \cdot y\_n \tag{7}$$

In a general case the values of *tstart* and *tend* can be different for different predictors. The above form of relationship between the predictor value and illumination is simple. More complex functions can be used in specific cases and produce a network with higher accuracy. However, if the relationship between the input and the illumination time is represented by a complex function, then one can argue that part of the computation is done, not by the network, but by the specific form of this relationship.

Coupling between oscillators, indicated by arrows in Figure 3, is achieved by reactions that extend the original Oregonator model. I assume that the coupling is of the activatory type and occurs via the exchange of reactor activators between oscillators [49]. Let *Ui* denote the activator of the *i*th node. The interactions between the nodes #*k* and #*j* appear as the results of reactions involving the activators *Uk* and *Uj*:

$$\mathcal{U}I\_{j} + \mathcal{B}\_{j,k} \to \mathcal{U}I\_{k} + \mathcal{C}\_{k} \tag{8}$$

$$\mathcal{U}I\_k + B\_{k,j} \to \mathcal{U}\_j + \mathcal{C}\_j \tag{9}$$

with the reaction rate constants *kB*,*j*,*<sup>k</sup>* and *kB*,*k*,*j*.

The changes in concentrations of *Uk* and *Uj* as the result of reactions (8) and (9) are described by:

$$\frac{du\_j}{dt} = -k\_{\text{B,j,k}} b\_{j,k} u\_j \tag{10}$$

$$\frac{du\_k}{dt} \quad = -k\_{B,k,j} b\_{k,j} u\_k \tag{11}$$

and the changes in concentration of *Uj* as the result of reaction (3) by:

$$\frac{d\mu\_j}{dt} = -k\_D d\_j \mu\_j \tag{12}$$

In Equations (10)–(12), *bj*,*k*, *bk*,*<sup>j</sup>* and *dj* denote concentrations of *Bj*,*k*, *Bk*,*<sup>j</sup>* and *Dj*, respectively. We assume that these concentrations are high with respect to the concentrations of the activators involved, are the same for all oscillators, and remain almost constant during the network evolution. Therefore, they are not included in the model of network evolution. Let us introduce symbols *α<sup>j</sup>* and *βj*,*k*, defined as: *α<sup>j</sup>* = *kDdj* and *βj*,*<sup>k</sup>* = *kB*,*j*,*kbj*,*k*. Keeping in mind that values of *α<sup>j</sup>* and *βj*,*<sup>k</sup>* can be modified by concentrations of *Dj* and *Bj*,*k*, we can treat them as free parameters that can be easily adjusted. The values of *α<sup>j</sup>* and *βj*,*<sup>k</sup>* are included in the optimization procedure. Let us also notice that all information on the network geometry is included in the values of *βj*,*k*, because for non-coupled nodes *βj*,*<sup>k</sup>* = *βk*,*<sup>j</sup>* = 0.

On the basis of the above assumptions, we can write equations describing the time evolution of the network:

$$\frac{du\_j}{dt} = \frac{1}{\varepsilon}(u\_j - u\_j^2 - (fv\_j + \phi\_j(t))\frac{u\_j - q}{u\_j + q}) - (a\_j + \sum\_{i=1,m} \beta\_{j,i})u\_j + \sum\_{i=1,m} \beta\_{i,j}u\_i \tag{13}$$

$$\frac{d\upsilon\_{\dot{j}}}{dt} = \upsilon\_{\dot{j}} - \upsilon\_{\dot{j}} \tag{14}$$

where i,j represent the *i*th and *j*th oscillator, respectively, and *m* is the number of oscillators in the network. This set of equations is solved numerically with a fourth-order RK algorithm [55] and time step of *dt* = 0.0005.

I assume that the output information can be extracted from observation of the network evolution during the time interval *Z* = [0, *tmax*]. This assumption is essential. In many previous studies on chemical computing, it was assumed that the system reaches a specific steady-state that represents the answer. In the present approach, the output information can be read from the time evolution of a selected oscillator observed in a finite time interval, and what happens later is irrelevant to computation.

Let us notice that Equation (5) has a physical meaning for any value of *tillum*(*j*). If *tillum*(*i*) < 0, then *φi*(*t*) is small and the oscillator #*i* is active during the whole observation interval. When *tillum*(*k*) > *tmax*, then *φk*(*t*) is large and the oscillator #*k* is inhibited within *Z* and does not oscillate.

In previous papers on classification with networks of interacting chemical oscillators, the number of activator maxima observed within the time interval *Z* − {0, *tmax*} on a selected node (say #*j* ) was considered as the output [47,49]. However, there is a question: What is the maximum of activator? I assume that *uj*(*t*) should have a maximum in a strict mathematical sense. This means that if *uj*(*t*) has a maximum at *t*0, then there exists *ν* > 0 such that *uj*(*t*) < *uj*(*t*0) for all *t* ∈ [*t*<sup>0</sup> − *ν*, *t*<sup>0</sup> + *ν*] ⊂ *Z*. Therefore, if *uj*(*t*) is growing at the end of *Z*, then *uj*(*tmax*) is not regarded as a maximum. Moreover, the value of *uj*(*t*0) should be larger than the activator concentration in the part of the oscillation cycle when the catalyst is in its reduced form. Here, I assume that the threshold value is 0.05. To illustrate maxima counting, let us consider the network with parameters given in the first row of Table 1. Figures 4 and 5 show the time evolution of activator (red curves) and inhibitor (blue curves) at all nodes for two selected points of the flag. The green line marks the threshold for the activator maximum (0.05). The results in Figure 4 illustrate the time evolution of activator and inhibitor on all nodes when the coordinates of the input point are: (−0.25, 0.25). The point characterized by such coordinates is located in the red area of the flag. For this input, there are two maxima of the activator at all network nodes. It can also be seen that on all nodes, the concentration of the activator is larger than 0.05 at the final stage of the evolution. However, the function *u*(*t*) is still growing at the end of *Z*, so the third maximum of *u*(*t*) is not observed. Figure 5 illustrates the time evolution of activator and inhibition on all nodes if the coordinates of the input point are: (0.29, −0.29). This point is located in the white region of the flag. Now, we observe three maxima of the activator at nodes #1 and #3, and two maxima at node #2. Therefore, if we consider the number of maxima on a selected oscillator as the network answer, then by selecting nodes #1 or #3, one can distinguish which of two points {(−0.25, 0.25),(0.29, −0.29)} is white and which is red. Such classification is not possible if node #2 is considered the output. Selection of the output node can be made on the basis of correct answers for the records included in the training dataset.

**Figure 4.** The time evolution of activator (the red curves) and inhibitor (the blue curves) observed on all nodes of the network defined by the parameters listed in the first line of Table 1. The coordinates of the input point are: (−0.25, 0.25). The green line marks the threshold for the activator maximum. There are 2 maxima at all nodes in the network.

**Figure 5.** The time evolution of activator (the red curves) and inhibitor (the blue curves) observed on all nodes of the network defined by the parameters listed in the first line of Table 1. The coordinates of the input point are: (−0.29, 0.29). The green line marks the threshold for the activator maximum. There are 3 maxima of *u*(*t*) on nodes #1 and #3 and two maxima of *u*(*t*) on node #2 within the observation time [0, *tmax*]. The green line marks the threshold for the activator maximum.

Figure 6 shows another example of classification with a number of activator maxima recorded on the selected node. Here, we consider the time evolution of the activator (the red curve) and the inhibitor (the blue curve) on node #1 of the network defined by the parameters listed in the second line of Table 1. Figure 6a illustrates the time evolution of concentrations for the point outside the red area (−0.25, 0.28), and Figure 6b for the point inside it (−0.39, −0.43). In the first case, no maximum of activator concentration is observed for times in *Z* − {0, *tmax*}; in the second case, we have a single maximum. Therefore, by observing the number of maxima on node #1, one can distinguish a red point from a white point, and a binary answer (0 or 1) is given. A classifying network with such properties can be regarded as highly optimized because the answer is obtained within a time similar to a single period of oscillation.

**Figure 6.** The time evolution of the activator (the red curve) and the inhibitor (the blue curve) on node #1 of the network defined by the parameters listed in the second line of Table 1: (**a**) *u*1(*t*) and *v*1(*t*) for the point inside the red area (−0.25, 0.28); (**b**) *u*1(*t*) and *v*1(*t*) for the point outside the red area (−0.39, −0.43).

Figures 7 and 8 illustrate another method for extracting the output information from the network evolution. This was inspired by [48], in which the output of an informationprocessing chemical BZ oscillator was related to the concentration of a selected reagent integrated over a specific time interval. Here we use the concentration of activator *u*(*t*) (red, Figure 7) and the concentration of inhibitor *v*(*t*) (blue, Figure 8) integrated over the observation time [0, *tmax*] as the network output. In these figures, the shaded areas below the function represent the integrals *Ju* = *tmax* <sup>0</sup> *<sup>u</sup>*1(*t*)*dt* and *Jv* <sup>=</sup> *tmax* <sup>0</sup> *v*3(*t*)*dt*, respectively. The integral value *J* is a real number, whereas an integer output is expected. In order to get such an output, I applied the following transformation:

$$output = floor(40 \cdot J)\tag{15}$$

where *J* stands for *Ju* or *Jv*. For the integrated concentration of the activator (Figure 7), we obtained *Ju* = 0.101 for the considered point inside the red area and *Ju* = 0.073 for the point outside it. Therefore, the network answers were 4 and 2 for these points, respectively. The method applied to the integral of the inhibitor (Figure 8) produced the outputs *Jv* = 0.203 for the considered point inside the red area and *Ju* = 0.181 for the point outside it. In this case, the network outputs were 8 and 7, respectively. Therefore, the integrals of the activator and inhibitor can also be used to classify points of different colors on the Japanese flag.

**Figure 7.** The time evolution of the activator at node #1 of the network defined by the parameters listed in the third line of Table 1: (**a**) *u*1(*t*) for the point (−0.25, 0.28) located inside the red area; (**b**) *u*1(*t*) for the point (−0.39, −0.43) located outside the red area. The red shaded area below the function represents the integral of *Ju* = *tmax* <sup>0</sup> *u*1(*t*)*dt*, considered as the network output.

**Figure 8.** The time evolution of the inhibitor at node #3 of the network defined by the parameters listed in the fourth line of Table 1: (**a**) *v*3(*t*) for the point (−0.25, 0.28) located inside the red area; (**b**) *v*3(*t*) for the point (−0.39, −0.43) located outside the red area. The blue shaded area below the function represents the integral *Jv* = *tmax* <sup>0</sup> *v*3(*t*)*dt*, considered as the network output.


**Table 1.** The parameters of networks that give the best correlations between the time evolution of the output oscillator and the point color.
