*2.4. Top-Down Design of Computing Networks*

To define an information-processing chemical oscillator network, we have to specify many parameters:


Obviously, a network with randomly selected parameters has a small chance of working as a good classifier. All parameters listed above should be optimized for executing the required function. I do not know any algorithm that allows for a straightforward design of the optimum oscillator network for a given problem. Still, we can apply the idea of supervised training for parameter optimization. Training means that we need a teacher, and in our optimization, it is a specific database *TA* that contains a sufficient number of records related to the considered problem *A* [56].

The network should include the output oscillator whose time evolution is transformed into the network answer. Here I consider two types of answers. One is the number of activator maxima observed in the time interval [0, *tmax*]. Another is the integral of the activator or inhibitor concentration. If the network parameters are known, the output oscillator can be located.

In order to find which oscillator should be used as the output, one can calculate the mutual information *I*(*S*;*Oj*) [57] between the discrete random variable *S* of record types in the training dataset *TA* (*S* = {*sn*, *n* = 1, *N*}) and the discrete random variable *Oj* representing the output of the *j*th oscillator in the network when the predictors of *n*th database record are used as the network input (*Oj* = {*oj*(*n*), *n* = 1, *N*}). The mutual information *I*(*S*;*Oj*) can be calculated as:

$$I(\mathcal{S}; O\_{\rangle}) = H(\mathcal{S}) + H(O\_{\rangle}) - H(\mathcal{S}, O\_{\rangle}) \tag{16}$$

where *H*() is the Shannon information entropy [58], and the discrete random variable (*S*,*Oj*) = {(*sn*, *oj*(*n*)), *n* = 1, *N*}. The oscillator #*i* for which the mutual information between *S* and *Oi* is maximal is used as the network output. The mutual information calculated for the output oscillator is considered the measure of network fitness:

$$Fitness = max\_{j \in \{1, 2, 3\}} I(S; O\_j) \tag{17}$$

It can be expected that in the majority of cases, optimization based on mutual information leads to a classifier with the highest accuracy [59].

The fitness function based on mutual information does not require specifying how the network evolution translates into the network output. However, if we know how to link the output *oj*(*n*) with the record type, then we can calculate the accuracy on node #*j* for training dataset *Ej* as the ratio between correctly classified cases and all cases of *TA*. For example, we can relate a specific number of activator maxima to the color of a point based on the majority of cases obtained for the training dataset. Therefore, we can alternatively define *Fitness* as:

$$Fitness = max\_{j \in \{1, 2, 3\}} E\_j \tag{18}$$

Both formulae (17) and (18) allow us to locate the output oscillator in the network.

It has been demonstrated that evolutionary optimization [60] oriented towards obtaining the best classifier for a representative training dataset of the problem can lead to a computing network that performs the anticipated task with reasonable accuracy [45,46,49]. In this approach, we represent the parameters as a code that undergoes recombination between the fittest members and mutations of an offspring. For recombination, two networks are selected, and their parameters are randomly separated into two parts. Next, an offspring is generated by combining one part of the first network with the other part of the second one. At this step, the function of an oscillator (input, normal) and illumination times of normal oscillators are copied to the offspring. For the next step, mutation of all parameters of the newborn offspring is considered. The probability of mutation rate is 0.5 per step, and the change in parameter value does not exceed 20%. Each generation of networks consists of the same number of elements. It includes a few fittest networks from the previous generation that are copied without changes in parameters. The remaining members of the next generation are offspring created by recombination and mutation operations applied to oscillators from the top 50% of networks from the previous population.

However, optimization of all parameters as mentioned above represents a computational problem of very high complexity. Before starting optimization, we introduce a number of simplifications:


After all these simplifications, the network is fully defined by the parameters listed in Table 1.

#### **3. What Is the Color of a Point on the Japanese Flag? (As Seen by the Networks)**

Here I present the results of network optimization. The training dataset is composed of 1000 records; 501 represent points inside the sun area and 499 are points outside it. The location of points can be seen, for example, in Figure 9. The applied evolutionary algorithm is a standard one, and it has been described in [45,46,49]. Optimization starts with 80 networks with randomly selected parameters. All networks in the initial population are evolved for maximum fitness. The fittest 10% (8) of the networks are copied to the next generation. The population of the next generation is completed to 80 networks by classifiers obtained by recombination and mutation operations applied to networks randomly selected from the upper 50% of the fittest. Next, the evolution step is repeated on the new population. The optimization procedure is executed for 1000 steps.

Figures 9–12 illustrate the answer of the fittest networks obtained for two Oregonator models and three different methods of interpreting the network output. The parameters of the discussed networks are given in Table 1. In all cases, the network output is based on the time evolution of a selected node in the network. For the first and second networks (Figures 9 and 10), this is the number of activator maxima observed on a selected node in the time interval [0, *tmax*]. For the third and fourth network, this is the integrated concentration of activator on a selected node (Figure 11) and the integrated concentration of inhibitor on a selected node (Figure 12). In all networks, the classification rule is derived after network optimization.

**Figure 9.** The answer of the network defined by the parameters listed in the first line of Table 1 to the records of the training dataset. Subfigures (**a**,**c**,**d**) are probability distributions of obtaining a given number of activator maxima on nodes #1, #2 and #3, respectively. The red bars correspond to points inside the red area; the blue bars refer to points outside the red area. Subfigure (**b**) illustrates correctly (yellow and red) and incorrectly (green and blue) classified points of the training dataset when node #1 is used as the output.

**Figure 10.** The answer of the network defined by the parameters listed in the second line of Table 1 to the records of the training dataset. Subfigures (**a**,**c**,**d**) are probability distributions of obtaining a given number of activator maxima on nodes #1, #2 and #3, respectively. The red bars correspond to points inside the red area; the blue bars refer to points outside the red area. Subfigure (**b**) illustrates correctly (yellow and red) and incorrectly (green and blue) classified points of the training dataset when node #1 is used as the output.

**Figure 11.** The answer of the network defined by the parameters listed in the third line of Table 1 to the records of training dataset. (**a**) The probability distribution of obtaining the value of *Ju* = *tmax* <sup>0</sup> *u*1(*t*)*dt* in the intervals [*k* ∗ 0.025,(*k* + 1) ∗ 0.025) for *k* ∈ {1, 2, 3, 4}. The red bars correspond to points inside the red area; the blue bars refer to points outside the red area. Subfigure (**b**) illustrates correctly (yellow and red) and incorrectly (green and blue) classified points of the training dataset.

Figure 9 shows the probability distribution of the number of activator maxima for the network defined by the parameters listed in the first line of Table 1. The red and blue bars correspond to points located in the sun and outside it, respectively. On nodes #2 and #3, we observe a large number of records representing points both in and out of the sun area that produced two maxima of the activator (cf. Figure 9c,d). On the other hand, for node #1, the number of activator maxima is two for most of the records representing points inside the sun area, and three maxima are observed for a large majority of points outside it. Therefore, node #1 is considered the output, and the classification rule is:


Such a rule incorrectly classifies 24 points located outside the sun and 8 points located inside the sun area. Therefore, the accuracy of this network with this classification rule on the training dataset is 96.8%. The geometry of correctly and incorrectly classified records of the training database is illustrated in Figure 9b. In this figure, as well as on similar figures for the other considered classifiers, the red dots represent points inside the sun area that are correctly recognized by the network. Similarly, the yellow dots represent points outside the sun area that are correctly recognized by the network. The green dots are points located within the sun area, but the network classifies them as belonging to the white region. Finally, the blue dots are misclassified points located in the white region. In order to see the errors more clearly, the dots corresponding to incorrectly classified points are larger than those representing correct answers. The location of errors at the edge of the red area is non-homogeneous and suggests the existence of regions where errors dominate.

**Figure 12.** The answer of the network defined by the parameters listed in the fourth line of Table 1 to the records of training dataset. (**a**) The probability distribution of obtaining the value of *Jv* = *tmax* <sup>0</sup> *v*1(*t*)*dt* in the intervals [*k* ∗ 0.025,(*k* + 1) ∗ 0.025) for *k* → {1, 2, 3, 4}. The red bars correspond to points inside the red area; the blue bars refer to points outside the red area. Subfigure (**b**) illustrates correctly (yellow and red) and incorrectly (green) classified points of the training dataset. As in (**a**,**b**) but for node #3: (**c**) the probability distribution of obtaining the value of *Ju* = *tmax* <sup>0</sup> *v*3(*t*)*dt* in the intervals [*k* ∗ 0.025,(*k* + 1) ∗ 0.025) for 3 ≤ *k* ≤ 10. Subfigure (**d**) illustrates correctly (yellow and red) and incorrectly (green and blue) classified points of the training dataset.

Similar results, but for the network defined by the parameters listed in the second line of Table 1, are shown in Figure 10. Again, node #1 is the network output because on nodes #2 and #3, we observe one activator maximum for a large number of records representing points located both in and out of the sun area (cf. Figure 10c,d). For node #1 (Figure 10a), most records representing points inside the sun area produced a single maximum, and for

most points located outside the sun, no maximum was observed (cf. Figure 6b). Therefore, node #1 is the output, and the classification rule is:


Such a rule incorrectly classifies just three points located outside the sun and seven points located inside the sun area. The accuracy of the second considered network with the classification rule listed above is 99% on the training dataset. The geometry of correctly and incorrectly classified records of the training database is illustrated in Figure 10b. All incorrectly classified points are located close to the boundary between the sun and the white region. Let us notice that the accuracy is 2% higher than for the previous network, and the only difference between these networks is the character of oscillations (cf. Figure 2b,c). This result suggests that "softer" oscillations may give networks with higher accuracy. This is the reason why the Oregonator Model II was used to investigate the alternative methods for relating the time evolution of concentration on a selected node with the network output.

Figure 11 presents results for the network defined by the parameters listed in the third line of Table 1. The network output is defined as: *Ju* = *tmax* <sup>0</sup> *uj*(*t*)*dt*, where *j* denotes the output node. The range of *Ju* values is divided into subintervals *Ik* with length 0.025 as follows: *Ik* = [*k* ∗ 0.025,(*k* + 1) ∗ 0.025) for *k* ∈ 1, 2, 3, 4. Figure 11a illustrates the probability distribution of *Ju* values for output node #1. The classification rule can be formulated as follows:


Such a rule incorrectly classifies just 2 points located inside the sun area and 11 points located outside it. All of the misclassified points are located close to the sun edge (cf. Figure 11b). The accuracy of the third of the considered networks with the above classification rule is 98.7% on the training database.

Keeping in mind the problem symmetry, we expect that node #1 should be the output. This was confirmed in the first three considered networks. However, node #1 gives poor accuracy for the network defined by the parameters listed in the fourth line of Table 1. For this network, the output is defined as: *Jv* = *tmax* <sup>0</sup> *vj*(*t*)*dt*, where *j* denotes the output node. As for the previously discussed network, the range of *Jv* values is divided into subintervals *Ik*. Figure 12a illustrates the probability distribution of *Jv* values for output node #1. The records representing points in the sun area and those located outside it return values of *Jv* in the third subinterval *I*3. Figure 12b shows the distribution of correctly and incorrectly classified records. For the majority rule derived from Figure 12a, all points located outside the sun are correctly classified, but there are many points inside the sun that are classified as located outside. Figure 12c illustrates the probability distribution of *Jv* values for output node #3. We observe a nice separation of records representing points in the sun area from those located outside. On the basis of this result, we formulate the following classification rule:


Such a rule incorrectly classifies just five points located inside the sun area and nine points located outside it. All of the misclassified points are located close to the sun's edge (cf. Figure 12d). The accuracy of the fourth considered network with the above classification rule is 98.6% on the training database. Network symmetry is reflected by high classification accuracy based on the inhibitor's time evolution on node #2. If we use this node as the output, we get classification accuracy of 97.9%. It can be expected that the difference between the accuracy with node #2 versus node #3 as the output is related to the choice of the training dataset, and in a perfectly balanced choice, both of these nodes should have the same accuracy.

#### **4. Discussion and Conclusions**

In the previous section, I demonstrated that three-node networks can be optimized to determine the color of a point in the Japanese flag with high accuracy. However, the accuracy was achieved on the training dataset, raising the question of whether the results shown in Figures 9–12 are general, or if they reflect correlations inevitable in the small set of randomly selected points that formed the training database.

To produce a stronger argument for the accuracy of optimized networks, I verify them on a large test dataset that contained 100,000 records. The points are generated randomly, and the test dataset includes 49,916 records corresponding to points inside the sun area and 50,084 points outside it.

The answers of the networks defined by the parameters listed in Table 1 to the records of the testing dataset are presented in Figure 13. Red and yellow dots represent correctly classified points located inside and outside the red area of the flag. To reduce the figure size, the number of points is limited to 10,000. The green color marks points located in the red area that are wrongly classified as being located in the white part. The blue color denotes points from outside that are classified as belonging to the sun area. In all cases, the classification is accurate, and the flag is nicely represented.

Figure 13a shows the results for the networks defined by the parameters listed in the first line of Table 1. As suggested by the results in Figure 9, the number of activator maxima on node #1 is used as the output. The majority rule is:


Such a rule incorrectly classifies 3216 points located outside the sun (blue in Figure 13a) and 1111 points located inside the sun area (green in Figure 13a). Therefore, the accuracy of this classifier is 95.668%. The sun area considered in this paper is more symmetrically oriented in the coordinate system than the Japanese flag located in the unit square [0, 1] × [0, 1] studied in [47]. For such a flag location the optimized classifier is characterized by: *tmax* = 20.23, *tstart* = 3.78, *tend* = 12.10 and *tillum*(3) = 6.37. Its accuracy tested on 100,000 records is 95.145%. However, to obtain this result, it was assumed that the values of parameter *α<sup>j</sup>* were identical for all nodes (= 0.849). Moreover, the values of *βi*,*<sup>j</sup>* = 0.251 were the same for all pairs of nodes. These assumptions definitely decreased the accuracy of the network reported in [47] if compared to the one defined in the first row of Table 1. Therefore, I expect that the accuracy of the three-oscillator classifier that codes the output in the number of activator maxima on a selected node does not systematically depend on the location of the flag in the coordinate system. The points classified as the sun (red and blue Figure 13a) form a characteristic horned disk already observed for a three-node classifier optimized to recognize the Japanese flag in the unit square [0, 1] × [0, 1] [47]. The fact that the horned shape is repeated by optimization using a training database containing points with coordinates in different ranges suggests that this shape is related to the parameters of the Oregonator Model I.

**Figure 13.** The answer of networks defined by the parameters listed in Table 1 to a testing dataset of 100,000 records. Yellow and red points are classified correctly. The network gives a wrong answer on points marked green (they belong to the sun but are classified as located outside it) and points marked blue (they are located outside the sun but are classified as belonging to the red area). Subfigures (**a**–**d**) correspond to networks with parameters listed in lines 1–4 of Table 1, respectively. The accuracy of these networks is (**a**) 0.957, (**b**) 0.984, (**c**) 0.976 and (**d**) 0.979, respectively.

The horned shape is not observed if one still uses the number of activator maxima as the output but considers another set of Oregonator parameters (Model II, (cf. Figure 13b)). A trace of undeveloped horns can be seen on the decreasing diagonal, but the sun disk is well represented. The classification rule (cf. Figure 10a) incorrectly predicts the color for 446 points located outside the sun (blue in Figure 13b) and 1192 points located inside the sun area (green in Figure 13b). The classifier accuracy is 98.362%. This result shows

that the change in the oscillator model can improve the accuracy and modify the geometry of incorrectly classified points.

Figure 13c,d represent classifiers where the integrals of activator and inhibitor observed on the output node were used as the output. The analysis of output related to activator concentration is shown in Figure 13c. The points inside the sun were those for which *Ju* ≥ 0.1. This rule incorrectly predicts the color for 2021 points located outside the sun (blue in Figure 13c) and only 401 points located inside the sun area (green in Figure 13c). The corresponding accuracy is 97.578%.

In the case of output related to inhibitor concentration (Figure 12c,d), the points inside the sun were those for which *Jv* ≥ 0.2. This rule incorrectly predicts the color for 1409 points located outside the sun (blue in Figure 13d) and only 691 points located inside the sun area (green in Figure 13d). The corresponding accuracy is 97.900%. The use of node #3 as the output is suggested due to its having the highest accuracy obtained on the training dataset. However, slightly higher accuracy (97.963%) is obtained if node #2 acts as the output. This indicates that even for such a simple problem as the Japanese flag one, the size of the training dataset is important to obtain the correct rule of classification.

The comparison between the accuracy of described classifiers shows that the method based on the integration of activator or inhibitor concentration is as good as the output represented by the number of the corresponding maxima. In the method of maxima counting, there was a magic parameter—the threshold value for acceptance of maximum. The method based on integration does not need to involve such additional parameters to obtain the result and seems to be a promising candidate for future investigation on the computing potential of networks composed of interacting chemical oscillators.

**Funding:** This research was funded by the Institute of Physical Chemistry as the support for the research activity of the Laboratory of Complex Systems and Chemical Processing of Information.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** All presented data are available from the author.

**Conflicts of Interest:** The author declares no conflict of interest.

#### **Abbreviation**

The following abbreviation is used in this manuscript: BZ Belousov–Zhabotinsky

#### **References**

