The Concilium of Information Processing Networks of Chemical Oscillators for Determining Drug Response in Patients With Multiple Myeloma

It can be expected that medical treatments in the future will be individually tailored for each patient. Here we present a step towards personally addressed drug therapy. We consider multiple myeloma treatment with drugs: bortezomib and dexamethasone. It has been observed that these drugs are effective for some patients and do not help others. We describe a network of chemical oscillators that can help to differentiate between non-responsive and responsive patients. In our numerical simulations, we consider a network of 3 interacting oscillators described with the Oregonator model. The input information is the gene expression value for one of 15 genes measured for patients with multiple myeloma. The single-gene networks optimized on a training set containing outcomes of 239 therapies, 169 using bortezomib and 70 using dexamethasone, show up to 71% accuracy in differentiating between non-responsive and responsive patients. If the results of single-gene networks are combined into the concilium with the majority voting strategy, then the accuracy of predicting the patient’s response to the therapy increases to ∼ 85%.

1 Introduction

The medical therapy of the future will be based on drugs individually selected for patient needs. Individual approach is necessary because it has been observed that one standard treatment does not work for all patients with the same type of cancer (Mulligan et al., 2007). An example is the therapy for patients with multiple myeloma. The drugs as bortezomib or dexamethasone have shown a good anti-myeloma effect, and they has been approved for treatment for patients with multiple myeloma (Field-Smith et al., 2006;¹). But the literature reports suggest that myeloma consists of many variants with different molecular pathologies (Hideshima et al., 2004; Mulligan et al., 2007). For certain subtypes, the drugs mentioned above can be effective, and for others, they may not. We can hope that using the gene expression profiling one can determine if certain drugs will be effective on a particular patient or not (Lesko and Woodcock, 2004). Considering the challenges one can face with genomic analysis for each patient (Zhan et al., 2006), we present a method that can be used to predict the drug effectiveness if the expression values of selected genes are known.

Our study is based on the results published in (Mulligan et al., 2007), its Supplementary Material (²) and the data of clinical tests available on the related web page (³). We considered this important medical problem to demonstrate the power of chemical computers operating with nonlinear chemical processes (Adamatzky et al., 2005). There are zillions of chemical computers operating around us because information processing in living organisms is based on chemical reactions. Animals and humans, using their nerve systems and brains, can control complex life processes, create models of the environment they function and even develop self-awareness. It demonstrates that Nature-made chemical computers can perform very complex computational tasks at low energy consumption. However, the development of human-made chemical computers has not shown as spectacular progress as semiconductor microprocessors (Waldrop, 2016). The binary information coding combined with robust logic gates, perfectly suited for electronic computers (Feynman et al., 2000), does not work for chemical computers because the lifetime of reagents is short. Therefore, the bottom-up approach (Feynman et al., 2000) does not help to design efficient chemical information processing devices. There are a few examples of the high computing potential of a chemical medium and its ability for parallel processing, such as the prairie-fire algorithm for labyrinth search (Steinbock et al., 1995) or image processing with a photosensitive variant of oscillatory Belousov-Zhabotinsky reaction (Kuhnert et al., 1989). However, the number of such human-proposed, efficient algorithms is limited. An alternative to the cleverness of a human researcher is the top-down machine design, where the computing medium is optimized to perform a selected task. The design of a chemical McCulloch-Pits neuron (McCulloch and Pitts, 1943) is difficult and requires high precision in setting the medium parameters (Gorecka and Gorecki, 2006). In this respect, networks of interacting chemical oscillators seem to be an interesting candidate for chemical computers. The recent results (Adamatzky et al., 2011; Holley et al., 2011; Szymanski et al., 2011; Adamatzky et al., 2012) demonstrated that networks of chemical oscillators can be easily assembled in experiments and studied for around a day (Muzika and Górecki, 2022). Such networks can be optimized to perform classification tasks (Gizynski et al., 2017) and process information with the best possible use of the chemical medium. The top-down design allows to design oscillator networks that perform functions for which a straightforward algorithm does not exist, as the determination of the cancer type on the basis of medical tests (Gizynski and Gorecki, 2017a).

In this paper, we describe another application of chemical oscillator networks for a medically oriented problem. We present a method that can help to determine the outcome of multiple myeloma therapy with bortezomib or dexamethasone drugs. The method is supposed to determine the response of a patient with a given gene expression profile to the therapy with the drugs mentioned above. The expression values of genes listed in Table 1 were considered. The functions of selected genes can be found in (⁴). For example, the gene RPS7 we consider in the discussed example of a single gene classifier encodes a ribosomal protein that belongs to the S7E family of ribosomal proteins. The gene CXCL5 encodes a protein that is a member of the CXC subfamily of chemokines, which recruit and activate leukocytes. This protein is proposed to bind the G-protein coupled receptor chemokine (C-X-C motif) receptor 2 to recruit neutrophils, promote angiogenesis, and remodel connective tissues. It is believed to play a role in cancer cell proliferation, migration, and invasion. The gene SERP1 is predicted to be involved in the endoplasmic reticulum unfolded protein response and protein glycosylation and act within several processes, including multicellular organism aging, positive regulation of organ growth, and positive regulation of peptide hormone secretion.

TABLE 1

TABLE 1. Genes considered for determining the drug response and their range.

In our study, we used the gene expression values provided by the authors of reference (Mulligan et al., 2007). Total RNA was isolated using Qiagen RNAeasy kit (⁵), and the expression values were measured by using the microarray technique (⁶). The information on the detailed procedure of myeloma cells enrichment, producing the expression of genes in myeloma cells, and quality control metrics can be found in Document 1 listed at the end of the web page (⁷). The records of the training dataset containing the gene expression values were generated from the GSM files describing the clinical results that can be found at (⁸). We included our training dataset in the Supplementary Material as the database file Table 1.xlsx. It contains the gene expression values data (columns A-O) together with the therapy result (Q column, 0 for nonresponsive and 1 for responsive case). Moreover, the database file includes information on the drug used (S column, PS341 for bortezomib and DEX for dexamethasone) and on the corresponding name of the GSM file with patient identification (U column). We considered information on 239 clinical tests. There were 169 tests with bortezomib, of which 84 were nonresponsive and 85 responsive, as well as 70 tests with dexamethasone, of which 42 were nonresponsive and 28 responsive.

Our approach is based on processing gene expression values with a network of chemical oscillators optimized for correlations between a single gene expression value and the result of therapy. The information processing network takes gene expression values as the inputs. The network functionality is optimized for maximum correlations between the result of therapy and the number of oscillations observed on the output oscillator, which is regarded as the network answer. The ideal network should simultaneously process information coming from many genes. However, it has been observed that in the case of many inputs, the evolutionary optimization of an information processing chemical network is long and ineffective because there are many local optima. In order to simplify the numerical simulations, we restricted our attention to networks made of three oscillators that process the expression data of a single gene. A typical optimized network returns ∼ 68% accuracy of prediction if for a given gene expression level, the therapy using bortezomib or dexamethasone is an effective treatment for multiple myeloma or not. We combined answers of single gene networks and made a concilium of networks based on the majority voting strategy. If such a strategy is applied to 15 considered genes, then we can predict the patient response to the drug therapy with an accuracy of ∼85% measured on the training dataset.

The paper is organized as follows: in Section 2 we define the procedure of input data normalization, introduce a 3-oscillator network and describe its optimization method. In the next section, we present optimization results and give the parameters of the best networks that correlate the individual gene expression values with the drug effectiveness. These networks are used for the concillum deciding if the drug therapy can be effective.

2 The Clinical Data and Their Classification

In the available dataset R = {r_k, k = 1, 239} we have 239 patient records and each record r_k contains information on the expression values of i = 15 genes (e_i,k, i = 1, 15) listed in Table 1. Moreover, Table 1 gives the maximum (Max_k(e_i,k)) and the minimum (Min_k(e_i,k)) gene expression values for each gene in R. The outcome of the therapy z_k ∈ {0, 1}, where 0 stands for nonresponsive therapy and 1 denotes responsive one, can be regarded as the record type. Therefore, a record r_k has the form of 16-tuple: r_k = (e_1,k, …, e_15,k, z_k). The problem of deciding if the therapy can be effective for a given patient reduces to the problem of finding an algorithm that gives the correct record type z provided that the predictor values {e_i, i = 1, 15} are known. It can be noticed that on this test group of patients, two trivial algorithms: 1) always use drugs and 2) do not use drugs because they will not help lead the correct therapy results in 47 and 53% cases respectively.

The values of a single gene expression weakly correlate with the efficiency of the drug therapy. Figure 1 illustrates the correlations between the expression of the RPS7 gene and the responsive/nonresponsive results observed in clinical tests. The whole range of gene expression values is divided into 10 subintervals (bins) of the same length, and their ranges are given in Table 2. It can be noticed that if the gene expression value of RPS7 is within bins no. 1,2,3 or 10, then the probability of successful therapy is higher than its failure. On the other hand, if the expression of the RPS7 gene is above 4,445 but below 11,065 then it is more likely that multiple myeloma will not respond to the drug. Using this rule, we can plan the therapy with an accuracy of 62.3%, which is much higher than that of the trivial algorithms mentioned above. Of 239 cases included in the dataset R, we obtained 70 correctly determined nonresponsive cases and 79 correctly determined responsive ones. We also observed 56 wrongly determined nonresponsive cases and 34 wrongly determined responsive ones. In the following, we show how this accuracy can be improved with simple classification networks based on interacting chemical oscillators.

FIGURE 1

FIGURE 1. The histogram showing correlation between the expression of the RPS7 gene and responsive (red)/nonresponsive (blue) results of drug therapy. The range of gene expression values corresponding to each bin are given in Table 2.

TABLE 2

TABLE 2. The range of gene expression values in the bins considered in the histogram shown in Figure 1.

It can be seen in Table 1 that the expression values significantly differ between genes. In order to unify the data for each gene we normalized the set of gene expression values e_i,k for all 239 patients using the formula:

$p_{i,k}=\frac{e_{i,k}-Mi{n}_k\left(e_{i,k}\right)}{Ma{x}_k\left(e_{i,k}\right)-Mi{n}_k\left(e_{i,k}\right)}(1)$

As an example the data {p_13,k, k = 1, 239} are illustrated in Figure 2. For clearer presentation of data we show points {(p_13,k, y_k), k = 1, 239} where y_k is a random number in the interval [0, 1]. The random y-coordinates Y = {y_k, k = 1, 239} were used to make the distribution of the training data easier to visualize, but of course, only the x-coordinate has a medical meaning. The same set of random y-coordinates Y was used to define points illustrated in Figures 3, 5, 7. The red and blue colors differentiate responsive and nonresponsive results of drug therapy, respectively. It can be seen that points corresponding to responsive and nonresponsive cases are not separated by the x-coordinate, related to the expression of gene no.13.

FIGURE 2

FIGURE 2. The distribution of normalized values of the RPS7 gene expression corresponding to responsive (red) and nonresponsive (blue) results of drug therapy respectively. The values of p_13,k are represented by the x-coordinate of marked points. The y-coordinate was randomly generated to differentiate points.

FIGURE 3

FIGURE 3. The distribution of normalized values of correctly and wrongly determined therapies using the majority rule based on the histogram of Figure 1. Red and blue points mark correctly determined responsive and nonresponsive records respectively. Black and green points mark wrongly determined responsive and nonresponsive records. The values of p_13,k are represented by the x-coordinate of marked points. The y-coordinate was randomly generated to differentiate points and it is the same as in Figure 2.

If we apply the rule following from the histogram shown in Figure 1 to the dataset of normalized data Q = {q_k, k = 1, 239} where q_k = (p_1,k, …, p_15,k, z_k) we obtain the distribution of correctly and wrongly determined results of drug therapy illustrated in Figure 3.

In the following, we consider the classification of the dataset Q using a network of interacting chemical oscillators. The idea of such computing was presented in the number of our papers (Gizynski and Gorecki, 2017a; Gizynski et al., 2017; Gorecki and Bose, 2020). We assume that there is a factor that controls oscillators and can inhibit oscillatory behavior. Such assumption is supported by the properties of Belousov Zabhotinsky(BZ) reaction (Belousov, 1959; Zhabotinsky, 1964) that has been widely used as a medium for chemical computing (Tóth and Showalter, 1995; Adamatzky et al., 2005; Yoshikawa et al., 2009; Gorecki et al., 2015; Dueñas-Díez and Pérez-Mercader, 2019; Proskurkin et al., 2020). The BZ-reaction is an oscillatory catalytic process (Epstein and Pojman, 1994). Among its reagents, we can distinguish HBrO₂ acting as the reaction activator and Br⁻ ions that are reaction inhibitors. It has been observed that for specific catalysts (for example, the ruthenium complex Ru(bpy)₃), the reaction is photosensitive (Kuhnert, 1986; Kuhnert et al., 1989). The illumination with a blue light generates Br⁻ ions that suppress oscillations (Kádár et al., 1997). The Oregonator model with two variables u and v representing concentrations of HBrO₂ and the oxidized form of the catalyst respectively is described by the equations (Rovinskii et al., 1984; Adamatzky et al., 2005):

$\frac{du}{dt}=\frac{1}{\epsilon }\left(u-u^2-\left(fv+\varphi \left(t\right)\right)\frac{u-q}{u+q}\right)(2)$

$\frac{dv}{dt}=u-v(3)$

The time evolution of a medium where BZ-reaction proceeds are determined by the values of parameters: f, q, and ɛ. The parameter ɛ sets up the ratio of time scale between variables u and v, q is the scaling constant, and f is the stoichiometric coefficient. The time-dependent function ϕ(t) is related to the medium illumination. The values of parameters f, q and ɛ can be selected such that for small ϕ(t) there are oscillations in u and v, but for a large ϕ the system converges to a stable stationary state (Kádár et al., 1997; Gorecki et al., 2014). Therefore, the value of ϕ can be interpreted as the light intensity in the Ru-catalyzed BZ-reaction. It can be used as an external factor to suppress oscillations or restore them. In the following, we used a modified Oregonator model to describe the time evolution of the considered computing oscillator networks.

There are a few arguments for the selection of the 2-variable Oregonator model. First, it provides a more realistic description of an oscillator network based on BZ-reaction than the oversimplified event-based-model used in early studies on computing networks of oscillators (Gizynski and Gorecki, 2017a; Gizynski et al., 2017). For example, the Oregonator model takes into account the effect of combined excitation of an oscillator by a few neighbors, which is missing in the event-based-model. Second, the model is still computationally simple, and it allows to perform a complex evolutionary optimization involving a huge number of evaluations of network evolution. Moreover, despite its simplicity, it provides a better than qualitative description of many phenomena related to BZ-reaction. It correctly describes the oscillation period as a function of reagent concentration and also can be used to model non-trivial phenomena like the migration of a spiral in an electric field (Sutthiopad et al., 2014) or reaction of a propagating pulse to time-dependent illumination (Tanaka et al., 2007). Of course, a model with a larger number of variables gives a more realistic description of BZ-reaction but, on the other hand, requires a more precise model of interactions between oscillators.

An example illustrating the idea of a considered computing oscillator network is illustrated in Figure 4A (Gorecki and Bose, 2020). The network is formed by three coupled oscillators marked by circles. We assume that the output information can be extracted from the observation of the network evolution during the time interval [0, t_max]. More precisely, the output information is coded in the number of activator maxima that are higher than a threshold value (in our study, this is 0.05) that are observed within the time interval [0, t_max] on a selected oscillator of the network. We consider time-dependent illumination ϕ_j(t) of the oscillator #j in the form:

${\varphi }_j\left(t\right)=0.1\cdot \left(\right.1.001+\tanh \left(-10\left(t-t_{\mathit{illum}}\left(j\right)\right)\right)(4)$

FIGURE 4

FIGURE 4. (A) An example illustrating the idea of a computing oscillator network. (B) The structure a network that produces a high correlation between the RPS7 gene expression and the result of drug therapy. The symbol In₁₃ marks input oscillators. The rightmost oscillator is a normal one, and it is also the output oscillator what is indicated with a double circle. (C) Correlations between the success of therapy and the number of activator maxima observed on the output oscillator.

This functional dependence is the same for every oscillator; however, the value of parameter t_illum(j) differs between oscillators. At the beginning of evolution, the inhibiting factor is high, which means that the oscillator is in a stationary state. For times t > t_illum(j) oscillations on the jth oscillator appear. At long times the value of ϕ_i(j) approaches 0.0001. For such illumination, Eqs 2, 3 produce oscillations characterized by the period of 8.2 time unit.

The use of illumination time (or, in general, the inhibition time for an oscillator) t_illum to influence oscillators is inspired by our experiments in which oscillations in individual BZ-droplets were controlled by blue LEDs (Gizynski and Gorecki, 2017b). In these experiments, we used just two illumination intensities: a low one for which the droplet was oscillating and a high one inhibiting oscillations. The transitions between the steady state and oscillations predicted by the 2-variable Oregonator model were in qualitative agreement with observations. For further studies, it will be interesting to consider more complex forms of time-dependent illumination because it has been observed that the time evolution of BZ-medium depends on the rate of changes in applied illumination (Tanaka et al., 2007).

Oscillators that form a computing network are of one of two types: the input oscillators and the normal ones (Gizynski and Gorecki, 2017a; Gizynski et al., 2017; Gorecki and Bose, 2020). If the jth oscillator is considered as a normal one, then the value of t_illum(j) is fixed. If this oscillator is consider as the input of a predictor p_i then the value of t_illum(j) is functionally related to p_i. For our analysis, we assume that the function has the form:

$t_{\mathit{illum}}\left(j\right)=t_{\mathit{start}}+\left(t_{\mathit{end}}-t_{\mathit{start}}\right)\ast p_i(5)$

The transformation given by Eq. 1 with parameters listed in Table 1 normalized the data included in the test dataset R to the interval [0, 1]. But we would like to know if the therapy is effective for any potential patient. Assume that the measured gene expression values are ${\tilde{e}}_i$, i = 1, 15. It can happen that some of ${\tilde{e}}_i$ are outside of the range listed in Table 1. Applying transformation (Mulligan et al., 2007) we obtain the corresponding ${\tilde{p}}_i$ outside [0, 1]. Nevertheless, according to Eq. 5 such input values also produce meaningful values of t_illum(j) that are accepted by gene information processing oscillator networks discussed below.

There is the output oscillator in the network. For fixed parameters describing the network, we select the output oscillator as the one that produces the highest accuracy of predicting record types for the dataset used for network training. In order to determine the network accuracy, we applied the following method. The network time evolution is simulated for all records of Q. For each oscillator and for each number of activator maxima, we formulate the relationship between the number of activator maxima and the outcome of therapy based on the majority of cases. The oscillator for which the number of errors is minimized is regarded as the output one.

The coupling between oscillators, indicated by two direction arrows in Figure 4A is achieved by reactions that extend the original Oregonator model. We assume that the coupling is of the activatory type and occurs via the exchange of reactor activators between oscillators (Gorecki and Bose, 2020; Bose and Gorecki, 2022). Let U_i denotes the activator of the ith oscillator. The exchange is described by reactions:

$U_j+B_j\to U_i+C_i(6)$

$U_i+B_i\to U_j+C_j(7)$

thus the coupling between oscillators is symmetric. We also assumed that the activator of each reaction can spontaneously decay in the process:

$U_j+D_j\to products(8)$

The symbols B, C and D appearing in the reactions above denote other molecules involved in these reactions. The concentrations of B and D (b, d) were assumed to be high with respect to the concentrations of activator and inhibitor, and hence their concentration was treated to be constant. If k_B and k_A are the reaction rate constants of reactions corresponding to coupling and decay respectively then the decrease in activator concentrations is described by the terms: k_Bb_ju_j, k_Bb_iu_i, and k_Ad_ju_j respectively. Having in mind high concentrations of B_j and D_j we can write those as βu_j, βu_j and αu_j where α and β are parameters with values controlled by b and d.

On the basis of the above assumptions we can formulate the following equations describing the time evolution of the network:

$\frac{d{u}_j}{dt}=\frac{1}{\epsilon }\left(u_j-u_j^2-\left(f{v}_j+{\varphi }_j\left(t\right)\right)\frac{u_j-q}{u_j+q}\right)-\left(\alpha +\beta \sum _{i=1,m}{s}_{j,i}\right)u_j+\beta \sum _{i=1,m}{s}_{j,i}{u}_i(9)$

$\frac{d{v}_j}{dt}=u_j-v_j(10)$

where i,j represent the jth and ith oscillator and m is the number of oscillators in the network. The variables u_j and v_j denote the concentration of an activator U_j and an inhibitor V_j respectively. The symbols s_j,i are defined as:

s_j,i = 0 if j = i or if j ≠ i and oscillators #j and #i do not interact,

s_j,i = 1 if j ≠ i and oscillators #j and #i do interact.

In our simulations we used the following set of Oregonator model parameters: ɛ = 0.3, q = 0.002, f = 1.1. For ϕ_j(t = 0) = 0.2 and these parameter values the stable steady state of Eqs 2, 3 is u_j = 0.00204 and v_j = 0.00204.

In order to define an information processing chemical oscillator network we have to specify many parameters: the number of oscillators m, the geometry of their connections (s_j,i), location of input and normal oscillators, all parameters for a model of chemical oscillations (ɛ, q and f), rates for reactions responsible for interactions between oscillators (α, β) the observation time t_max, illumination times for all normal oscillators t_illum(i) and parameters t_start and t_end that translate an input value into the illumination of an input oscillator (cf. Eq. 5). The problem is even more complicated as we can consider other models of chemical oscillators than Oregonator, different functions linking input values with the time interval within which the inhibiting factor is applied, and various models for coupling between oscillators. We do not know any algorithm that allows for a straightforward design of the optimum oscillator network for a given problem. Still, we can apply a parameter optimization algorithm to a training dataset in the hope it produces a network that gives a reasonable solution to the problem. However, optimization of all parameters mentioned above represents a computational problem of very high complexity. Before starting the optimization, we introduced a number of simplifications:

(1) we restricted our attention to classifiers formed by m = 3 oscillators,

(2) we assumed that each oscillator interacted with all others, so s_j,i ≡ 1 for j ≠ i. The geometry of such network is illustrated in Figure 4A.

(3) There has to be an input oscillator in the network and a normal one. Without the input oscillator, the network returns the same answer on all inputs. Without the normal oscillator, the network evolves like a single oscillator. Keeping in mind the symmetry of the considered network, we can assume that the oscillator #1 is the input oscillator and the oscillator #3 is a normal one. The role of the oscillator #2 is the subject of optimization.

After these simplifications the network is fully characterized by t_max, α, β, t_start, t_end, t_illum (3), the role of oscillator #2 and if it is the normal one, its illumination time t_illum (2). We optimized the values of these parameters using Q as the training dataset. The fitness function of evolutionary optimization was the maximum accuracy between the network output coded in the number of activator maxima observed on one of the oscillators and the record type z (Gizynski and Gorecki, 2017a; Gizynski et al., 2017; Gorecki and Bose, 2020). The time evolution of networks was obtained by solving Eqs 9, 10 numerically using 5-th order Cash-Karp algorithm (Cash and Karp, 1990) with dt = 10^–3 time steps. For a given network, all three oscillators were considered as potential candidates for output one. The oscillator that produced the highest accuracy on the training dataset was regarded as the output one. The applied evolutionary algorithm is a standard one (Goldberg, 1989), and it has been described in our previous papers (Gizynski and Gorecki, 2017a; Gizynski et al., 2017; Gorecki and Bose, 2020). The optimization started with 100 networks with randomly generated parameters. The next generation of networks included 5 top fit networks of the previous generation and 95 networks formed by recombination of parameters of two networks selected from 40 best networks of the previous generation. Each network obtained by recombination was then allowed to mutate. Mutations included the values of all parameters and the type of oscillator #2. The optimized network was obtained after 500 evolutionary steps.

3 Results

The evolutionary optimization described in the previous Section was used to design networks with a high correlation between the number of activator maxima observed on one of the oscillators and the result of drug therapy. As an example, we show such a network for linking the RPS7 normalized gene expression value with the success of therapy. For this input, the evolutionary optimization produced the network illustrated in Figure 4B. It is composed of two input oscillators, marked with In₁₃, that accept the value of p_13,k. The rightmost oscillator, marked with the double circle, is the output one. The circle marking this oscillator is also a base for a pie chart representing the ratio t_illum (3)/t_max by the surface of the red slice. Figure 4C shows the distribution of a number of activator maxima observed on the output oscillator for responsive and nonresponsive multiple myeloma treatment. On this basis, we can define the rule that ensures the highest accuracy on the training dataset Q. For this network, the rule is:

- if 1,3,4,5,9 or 10 activator maxima are observed on the output oscillator, then the patient belongs to the responsive group.

- if another number of activator maxima is observed, then unsuccessful treatment with the drugs is expected.

If applied to the training dataset Q, this rule gives 71.1% of correct answers, which is 5% higher than the trivial rule based on the distribution of gene expression values (cf. Figure 1). Of 239 cases included in the dataset Q, we obtained 91 correctly determined nonresponsive cases and 79 correctly determined responsive ones. We also observed 35 wrongly determined nonresponsive cases and 34 wrongly determined responsive ones. The distribution of correctly and incorrectly classified points from the training dataset is illustrated in Figure 5. Here again, the y-coordinate is random, and it is the same as introduced to differentiate the values of p_13,k in Figure 2.

FIGURE 5

FIGURE 5. Location of correct and incorrect predictions on the result of drug therapy by a network of oscillators optimized for the RPS7 gene. The values of p_13,k are represented by the x-coordinate of marked points. The y-coordinate is randomly generated to differentiate points and is the same as in Figure 2.

We optimized 3-oscillator networks for all genes listed in Table 1. The parameters of all optimized networks are listed in Table 3. The gene numbers in Table 3, Table 4, and Figure 6 correspond to those in Table 1. The rules that translate the number of activator maxima on the output oscillator and the therapy result are given in Table 4. The accuracy of optimized information processing networks is shown in Figure 6, and in Table 4. It is in the range between 66.5% (genes SERP1, CXCL5 and IL15) to 71.1% (gene RPS7).

TABLE 3

TABLE 3. The parameters of networks that give the best correlations between the number of activator maxima on the output oscillator and the therapy result.

TABLE 4

TABLE 4. The rules that translate the number of activator maxima on the output oscillator and the effective therapy using bortezomib or dexamethasone drugs.

FIGURE 6

FIGURE 6. Comparison between the accuracy of predictions of the result of drug therapy based on the histogram of gene expression values (red bars) and the optimized network of oscillators (blue bars). The gene numbers correspond to these in Table 1.

To increase the accuracy in determining the success or failure of drug therapy, we called a concilium of optimized networks. Each member of concilium is a network specialized in finding correlations between the expression value of one gene and the success of drug therapy and has one vote. The final decision is taken on the basis of majority voting. The accuracy of determining if the therapy using the bortezomib or dexamethasone drugs is efficient or not for different majority rules applied to the concilium is shown in the second column of Table 5. We can see that the decision based on the opinions of more than half of concilium members is accurate for almost 85% of cases included in the training dataset. Of 239 cases included in the dataset Q we obtained 117 correctly determined nonresponsive cases and 86 correctly determined responsive ones. We also observe 9 wrongly determined nonresponsive cases and 27 wrongly determined responsive ones. Their distribution is illustrated in Figure 7.

TABLE 5

TABLE 5. The accuracy of determination if the therapy using the bortezomib or dexamethasone drugs is efficient for different majority rules. The results for the optimized network and the network with modified parameters (cf. Table 6) are compared.

FIGURE 7

FIGURE 7. Location of correct and incorrect predictions on the result of drug therapy predicted by a concilium of networks of oscillators optimized for all genes. The values of p_13,k are represented by the x-coordinate of marked points. The y-coordinate is the same as in Figure 2.

But do we need oscillator networks? Correlations between the values of single gene expression corresponding to responsive and nonresponsive cases from the training dataset and the results of drug therapy can be extracted from the histograms of single gene expression values, as it was done for the RPS7 gene (cf. Figure 1). For this gene, the accuracy of such a method is 62.3% which is not much lower than that of the optimized classifier (71.1%). So why not try to make the conciluim based on histograms for all genes? We have tested such an approach by dividing the whole range of gene expression values into 10 subintervals and introducing the rule based on the majority of cases in subintervals for each gene. Next, we applied the majority voting strategy for all records of Q. The accuracy of such concilium was 66.9%. Of 239 cases included in the dataset Q, we obtained 100 correctly determined nonresponsive cases and 60 correctly determined responsive ones. We also observed 26 wrongly determined nonresponsive cases and 53 wrongly determined responsive ones. Therefore the accuracy of such concilium is much smaller than the concilium based on networks of oscillators optimized for correlations between the single gene expression value and the outcome of therapy.

4 Conclusion and Discussion

In this paper we discussed the application of information processing network formed by chemical oscillators for determination of the outcome of the multiple myeloma therapy with bortezomib or dexamethasone drugs. The network input information comes form the gene expression values. Information was processed by simple networks, each made of 3 oscillators. Each network was optimized to find correlations between the expression value of a particular gene and the outcome of the therapy. Individual classifiers gave the accuracy in the range between 66.5 and 71.1% (cf. Figure 7). To improve the determination of the therapy outcome we considered the concilium of 15 classifiers and accepted the majority decision. Such strategy increased the accuracy to almost 85%, which seems to be a promising result for the further development of the method.

In the mid columns of Table 5 we presented the accuracy of the concilium method based on classifiers optimized for the whole training dataset to therapies in which one of the drugs (bortezomib or dexamethasone) was used. The idea was to find if the therapy prediction accuracy depends on the drug. The results are similar, so we can conclude that for both drugs, the patient genetic profile is similarly correlated with the therapy success. It would be interesting to make a similar concilium separately for each drug, but the solution to this problem requires a much larger database of clinical trials than the one we had access to.

To check how sensitive to fluctuations of parameters are the results produced by the concilium formed of optimized networks, we considered random modifications in parameter values. For each optimized network, we selected one parameter at random and decreased or increased its value by ± 1%. The details on applied modifications and their influence on the accuracy of each network are given in Table 6. In all cases, the accuracy decreased by 2–3%. A similar decrease of accuracy is observed for the decision of concilium (82.8% for the concilium of modified networks). Still, such accuracy is high enough to claim that concilium strategy for determination of drug effectiveness is robust to random changes in parameters and fluctuations in the medium.

TABLE 6

TABLE 6. The accuracy of modified oscillator networks for correlations between the gene expression value and the therapy result. The table also defines modification introduced to the optimized network.

We can suggest two ways in which the accuracy in the determination of therapy effectiveness can be increased:

One of them is to consider the voting strategy with more complex networks formed by a larger number of oscillators that are used to determine correlations between a single gene expression value and the therapy outcome. One can expect that “wiser” members of concilium can produce more accurate answers. However, the strategy of employing top specialists does not guarantee top results. Our simulations have shown that the synergy between concilium members is also important and should be taken into account. We continued optimization for some networks processing gene expression values and obtained higher accuracy than that listed in Table 4. However, if we replaced the optimized network for the gene SERP1 that led to the accuracy of 66.5% (cf. #1 in Table 3) by a wiser member (accuracy 67.7%, t_max = 80, t_start = 3.52, t_end = 32.23, = 0.71, β = 0.09, normal oscillators 2 and 3, $t_{\mathit{illum}}^2=3.38$, $t_{\mathit{illum}}^3=4.24$) than the accuracy of concilium decreased to 83.6%.

Alternatively, one can consider networks that are processing expression values for more than a single gene. However, the pairs of gene expression values corresponding to responsive and nonresponsive therapies do not show clear separation in the square [0, 1] × [0, 1] (cf. Figure 8). Therefore, it can be anticipated that a large oscillator network is necessary for the data classification, and its optimization will be numerically complex.

FIGURE 8

FIGURE 8. The distribution of normalized expression values of the RPS7 and CFLAR genes corresponding to responsive (red) and nonresponsive (blue) results of drug therapy respectively. The point coordinates are (p_13,k, p_8,k), k = 1, 239.

Finally, let us make more general comments on the importance of the presented results:

First, they demonstrate a high computing potential of networks composed of interacting oscillators. The distributions of gene-expression values corresponding to responsive and non-responsive patients do not show a clear separation, as illustrated in Figure 2 and indicated by the histogram in Figure 1. We can expect that classification of these data with classical neural networks is inefficient and requires a large number of nodes. We demonstrated that the cases separation with a reasonable accuracy exceeding 65% can be done with a network of just 3 oscillators. Suppose the effectiveness of oscillator networks is confirmed on other problems. Then, as a natural development of this approach, we can expect a new class of integrated circuits made with semiconductors with easy control of the geometry of interactions and parameters of oscillators. They can operate similar to Intel Neural Compute Stick designed to support computation with classical neural networks⁹.

Second, the semiconductor devices work properly in a narrow range of temperatures around the room ones, whereas the range of conditions in which chemical oscillations are observed is much wider. Therefore, we can think of designing chemical computers for the specific environment they are supposed to function, for example, for space research applications. Furthermore, chemical computers operate on the energy of their reagents, so they do not need additional energy supply.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author Contributions

All authors listed have made a substantial, direct, and intellectual contribution to the work and approved it for publication.

Funding

This publication is part of a project that has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 711859. Scientific work funded from the financial resources for science in the years 2017–2022 awarded by the Polish Ministry of Science and Higher Education for the implementation of an international co-financed project.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fchem.2022.901918/full#supplementary-material

Footnotes

¹https://www.cancertherapyadvisor.com/home/cancer-topics/hematologic-cancers/hematologic-cancers-treatment-regimens/multiple-myeloma-treatment-regimens/

²https://ashpublications.org/blood/article/109/8/3177/23711/Gene-expression-profiling-and-correlation-with

³https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE9782

⁴https://www.ncbi.nlm.nih.gov/gene

⁵https://www.qiagen.com/us/products/discovery-and-translational-research/dna-rna-purification/rna-purification/total-rna/rneasy-kits/

⁶https://en.wikipedia.org/wiki/DNA_microarray

⁷https://ashpublications.org/blood/article/109/8/3177/23711/Gene-expression-profiling-and-correlation-with

⁸https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE9782

⁹https://www.intel.com/content/www/us/en/developer/tools/neural-compute-stick/overview.html

References

Adamatzky, A., De Lacy Costello, B., and Asai, T. (2005). Reaction–diffusion Computers. New York, NY, USA: Elsevier.

Google Scholar

Adamatzky, A., Holley, J., Bull, L., and De Lacy Costello, B. (2011). On Computing in Fine-Grained Compartmentalised Belousov-Zhabotinsky Medium. Chaos, Solit. Fractals 44, 779–790. doi:10.1016/j.chaos.2011.03.010

CrossRef Full Text | Google Scholar

Adamatzky, A., Holley, J., Dittrich, P., Gorecki, J., De Lacy Costello, B., Zauner, K. P., et al. (2012). On Architectures of Circuits Implemented in Simulated Belousov-Zhabotinsky Droplets. Biosystems 109, 72–77. doi:10.1016/j.biosystems.2011.12.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Belousov, B. P. (1959). Collection of Short Papers on Radiation Medicine. Moscow: Medgiz, 145–152.

Google Scholar

Bose, A., and Gorecki, J. (2022). Computing with Networks of Chemical Oscillators and its Application for Schizophrenia Diagnosis. Front. Chem. 10. doi:10.3389/FCHEM.2022.848685

CrossRef Full Text | Google Scholar

Cash, J. R., and Karp, A. H. (1990). A Variable Order Runge-Kutta Method for Initial Value Problems with Rapidly Varying Right-Hand Sides. ACM Trans. Math. Softw. 16, 201–222. doi:10.1145/79505.79507

CrossRef Full Text | Google Scholar

Dueñas-Díez, M., and Pérez-Mercader, J. (2019). How Chemistry Computes: Language Recognition by Non-Biochemical Chemical Automata. From Finite Automata to Turing Machines. iScience 19, 514–526. doi:10.1016/j.isci.2019.08.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Epstein, I. R., and Pojman, J. A. (1994). Introduction to Nonlinear Chemical Dynamics: Oscillations, Waves, Patterns, and Chaos. New York, NY, USA: Oxford University Press.

Google Scholar

Feynman, R. P., Hey, T., and Allen, R. (2000). Feynman Lectures on Computation. Boulder, Colorado, USA: CRC Press.

Google Scholar

Field-Smith, A., Morgan, G. J., and Davies, F. E. (2006). Bortezomib (Velcade?) in the treatment of multiple myeloma. Ther. Clin. Risk Manag. 2 (3), 271–279. doi:10.2147/tcrm.2006.2.3.271

PubMed Abstract | CrossRef Full Text | Google Scholar

Gizynski, K., and Gorecki, J. (2017). Cancer classification with a network of chemical oscillators. Phys. Chem. Chem. Phys. 19, 28808–28819. doi:10.1039/c7cp05655a

PubMed Abstract | CrossRef Full Text | Google Scholar

Gizynski, K., and Gorecki, J. (2017). Chemical memory with states coded in light controlled oscillations of interacting Belousov-Zhabotinsky droplets. Phys. Chem. Chem. Phys. 19 (9), 6519–6531. doi:10.1039/c6cp07492h

PubMed Abstract | CrossRef Full Text | Google Scholar

Gizynski, K., Gruenert, G., Dittrich, P., and Gorecki, J. (2017). Evolutionary Design of Classifiers Made of Droplets Containing a Nonlinear Chemical Medium. Evol. Comput. 25, 643–671. doi:10.1162/evco_a_00197

PubMed Abstract | CrossRef Full Text | Google Scholar

Goldberg, D. E. (1989). Genetic Algorithms in Search, Optimization and Machine Learning. Boston, MA: Addison-Wesley Longman Publishing Co., Inc.

Google Scholar

Gorecka, J., and Gorecki, J. (2006). Multiargument logical operations performed with excitable chemical medium. J. Chem. Phys. 124, 084101. doi:10.1063/1.2170076

PubMed Abstract | CrossRef Full Text | Google Scholar

Gorecki, J., and Bose, A. (2020). How Does a Simple Network of Chemical Oscillators See the Japanese Flag? Front. Chem. 8, 580703. doi:10.3389/fchem.2020.580703

PubMed Abstract | CrossRef Full Text | Google Scholar

Gorecki, J., Gizynski, K., Guzowski, J., Gorecka, J. N., Garstecki, P., Gruenert, G., et al. (2015). Chemical computing with reaction-diffusion processes. Phil. Trans. R. Soc. A 373, 20140219. doi:10.1098/rsta.2014.0219

PubMed Abstract | CrossRef Full Text | Google Scholar

Gorecki, J., Gorecka, J. N., and Adamatzky, A. (2014). Information coding with frequency of oscillations in Belousov-Zhabotinsky encapsulated disks. Phys. Rev. E 89, 042910. doi:10.1103/PhysRevE.89.042910

PubMed Abstract | CrossRef Full Text | Google Scholar

Hideshima, T., Bergsagel, P. L., Kuehl, W. M., and Anderson, K. C. (2004). Advances in biology of multiple myeloma: clinical applications. Blood 104, 607–618. doi:10.1182/blood-2004-01-0037

PubMed Abstract | CrossRef Full Text | Google Scholar

Holley, J., Adamatzky, A., Bull, L., De Lacy Costello, B., and Jahan, I. (2011). Computational modalities of Belousov-Zhabotinsky encapsulated vesicles. Nano Commun. Netw. 2, 50–61. doi:10.1016/j.nancom.2011.02.002

CrossRef Full Text | Google Scholar

Kádár, S., Amemiya, T., and Showalter, K. (1997). Reaction Mechanism for Light Sensitivity of the Ru(bpy)32+-Catalyzed Belousov−Zhabotinsky Reaction. J. Phys. Chem. A 101, 8200–8206. doi:10.1021/jp971937y

CrossRef Full Text | Google Scholar

Kuhnert, L. (1986). A new optical photochemical memory device in a light-sensitive chemical active medium. Nature 319, 393–394. doi:10.1038/319393a0

CrossRef Full Text | Google Scholar

Kuhnert, L., Agladze, K. I., and Krinsky, V. I. (1989). Image processing using light-sensitive chemical waves. Nature 337, 244–247. doi:10.1038/337244a0

CrossRef Full Text | Google Scholar

Lesko, L. J., and Woodcock, J. (2004). Translation of pharmacogenomics and pharmacogenetics: a regulatory perspective. Nat. Rev. Drug Discov. 3, 763–769. doi:10.1038/nrd1499

PubMed Abstract | CrossRef Full Text | Google Scholar

McCulloch, W. S., and Pitts, W. (1943). A logical calculus of the ideas immanent in nervous activity. Bull. Math. Biophysics 5, 115–133. doi:10.1007/BF02478259

CrossRef Full Text | Google Scholar

Mulligan, G., Mitsiades, C., Bryant, B., Zhan, F., Chng, W. J., Roels, S., et al. (2007). Gene expression profiling and correlation with outcome in clinical trials of the proteasome inhibitor bortezomib. Blood 109 (8), 3177–3188. doi:10.1182/blood-2006-09-044974

PubMed Abstract | CrossRef Full Text | Google Scholar

Muzika, F., and Górecki, J. (2022). Identification of the best medium for experiments on chemical computation with Belousov-Zhabotinsky reaction and ferroin-loaded Dowex beads. Reac Kinet. Mech. Cat. 135, 1187–1209. doi:10.1007/s11144-022-02171-4

CrossRef Full Text | Google Scholar

Proskurkin, I. S., Smelov, P. S., and Vanag, V. K. (2020). Experimental verification of an opto-chemical "neurocomputer". Phys. Chem. Chem. Phys. 22 (34), 19359–19367. doi:10.1039/d0cp01858a

PubMed Abstract | CrossRef Full Text | Google Scholar

Rovinskii, A. B., Zhabotinskii, A. M., and Irving, R. (1984). Mechanism and mathematical model of the oscillating bromate-ferroin-bromomalonic acid reaction. J. Phys. Chem. 88, 6081–6084. doi:10.1021/j150669a001

CrossRef Full Text | Google Scholar

Steinbock, O., Tóth, Á., and Showalter, K. (1995). Navigating Complex Labyrinths: Optimal Paths from Chemical Waves. Science 267, 868–871. doi:10.1126/science.267.5199.868

PubMed Abstract | CrossRef Full Text | Google Scholar

Sutthiopad, M., Luengviriya, J., Porjai, P., Tomapatanaget, B., Müller, S. C., and Luengviriya, C. (2014). Unpinning of spiral waves by electrical forcing in excitable chemical media. Phys. Rev. E 89 (5), 052902. doi:10.1103/PhysRevE.89.052902

PubMed Abstract | CrossRef Full Text | Google Scholar

Szymanski, J., Gorecka, J. N., Igarashi, Y., Gizynski, K., Gorecki, J., Zauner, K. P., et al. (2011). Droplets with information processing ability. Int. J. Unconv. Comput. 7, 185–200.

Google Scholar

Tanaka, M., Nagahara, H., Kitahata, H., Krinsky, V., Agladze, K., and Yoshikawa, K. (2007). Survival versus collapse: abrupt drop of excitability kills the traveling pulse, while gradual change results in adaptation. Phys. Rev. E 76 (1 Pt 2), 016205. doi:10.1103/PhysRevE.76.016205

PubMed Abstract | CrossRef Full Text | Google Scholar

Tóth, Á., and Showalter, K. (1995). Logic gates in excitable media. J. Chem. Phys. 103, 2058–2066. doi:10.1063/1.469732

CrossRef Full Text | Google Scholar

Waldrop, M. M. (2016). The chips are down for Moore's law. Nature 530 (7589), 144–147. doi:10.1038/530144a

PubMed Abstract | CrossRef Full Text | Google Scholar

Yoshikawa, K., Motoike, T. I., Yamaguchi, T., Igarashi, Y., Gorecki, J., and Gorecka, J. N. (2009). Basic information processing operations with pulses of excitation in a reaction-diffusion system. Int. J. Unconv. Comput. 5, 3–37.

Google Scholar

Zhabotinsky, A. M. (1964). Periodic liquid phase reactions. Proc. Acad. Sci. USSR 157, 392–395.

Google Scholar

Zhan, F., Huang, Y., Colla, S., Stewart, J. P., Hanamura, I., Gupta, S., et al. (2006). The molecular classification of multiple myeloma. Blood 108, 2020–2028. doi:10.1182/blood-2005-11-013458

PubMed Abstract | CrossRef Full Text | Google Scholar