Open access peer-reviewed chapter - ONLINE FIRST

Model of an Artificial Blastula for Assessing Development Toxicity

Written By

František Muzika and Jerzy Górecki

Submitted: 20 January 2023 Reviewed: 23 January 2023 Published: 04 March 2023

DOI: 10.5772/intechopen.110260

Bioinformatics and Medical Informatics - Annual Volume 2024 IntechOpen
Bioinformatics and Medical Informatics - Annual Volume 2024 Authored by Slawomir Wilczynski

From the Annual Volume

Bioinformatics and Medical Informatics - Annual Volume 2024 [Working Title]

Associate Prof. Slawomir Wilczynski

Chapter metrics overview

78 Chapter Downloads

View Full Metrics

Abstract

We are concerned with computer simulations of a ring of 20 coupled CSTRs with glycolytic oscillatory reaction. Each CSTR represents an artificial cell, and the ring can be regarded as an artificial blastula. The cells are coupled to two adjacent CSTRs via the mass exchange of reagents. The glycolytic oscillatory reaction is simulated using the two-variable core model. Our work is focused on the classification of stationary discrete nonuniform concentration patterns (discrete Turing patterns). The control parameters in simulations are autocatalytic and inhibition rate coefficients, as well as the transport rate coefficients. We performed the analysis of stability and bifurcations of stationary states to identify the stationary states. The inflow of reagents into each CSTR was used to initiate a particular pattern. We propose a method to assess the morphogenetic toxicity of any chemical from a database by switching between patterns or between patterns and oscillations. Moreover, we investigated nonuniform patterns that create discrete concentration waves inside the ring of 20 coupled cells, which can trigger gastrulation.

Keywords

  • discrete turing patterns
  • glycolysis
  • 2D artificial blastula
  • cheminformatics
  • morphogenesis toxicity

1. Introduction

Chemoinformatics tools allow us to define information and parameters of various chemicals influencing more complex systems. The information and parameters are focused on chemical structure, experimental and physicochemical data, and toxicity. The toxicity data can concern living tissues, metabolic pathways [1], DNA [2, 3], RNA, [4, 5], gut microbia [6], or mitochondrial toxicity [7]. Chemoinformatics is also essential for the development of new drugs [8, 9, 10, 11] and for the assessment of their toxicity [12]. Since the databases of chemicals are extensive and count in millions of elements, the approach of using artificial intelligence [13], artificial neural networks [14], and machine learning seems to be the only method to deal with such a huge amount of collected information [15].

In this paper, we are concerned with the toxicity of chemicals causing morphogenetic malfunctions. The malfunctions can be seen as deformed organs, excessive or missing limbs, and fingers. The toxic factors can also terminate morphogenesis, causing the death of living organisms. On the other hand, the identification of toxic factors can be useful for slowing down the growth of cancer cell clusters. As an alternative to experiments with actual living cells, we created the 2D model of an artificial blastula, which mimics a multicellular living organism in its developmental stage (blastula). Within our model, toxicity and morphogenic malfunctions can be recorded as the formation and destruction of discrete nonuniform patterns and switching between them. These patterns were first introduced by A. Turing in his pioneering study [16]. They are partially responsible for yeast budding [17] and are expected to play a key role during fetal development and gastrulation [18], limb development [19], and fingers development [20]. While the stationary Turing pattern can occur under conditions of faster transport rate of inhibitor species, usually under the scheme “short-range activation and long-range inhibition” [21, 22], they can also occur for equal transport rates of both activator and inhibitor, provided that the overall kinetics is enhanced, and the system is properly perturbed [23, 24, 25, 26]. Turing patterns have been proven experimentally [27, 28, 29, 30, 31, 32], but for its occurrence, the activator molecule transport rate has to be slowed down by different agents. The occurrence of discrete Turing patterns has been shown experimentally and in models in various artificial systems using the Belousov-Zhabotinsky reaction and its modifications [33, 34, 35, 36, 37, 38].

Our recent research shows it is possible to make transitions between uniform oscillations and discrete nonuniform patterns [39]. This can be used for chemical computing as well as a type of chemical memory. Moreover, nonuniform patterns can also be used as the output of chemical classification. The network considered in this paper can be used as a classifier under different kinetic and transport parameters and under different network topology. The reaction parameters and coupling constants have to be optimized for each topology. Evolutionary algorithms are frequently used for parameter optimization. Once optimized, even a small network can classify schizophrenia [40] or color of points on the Japanese flag [41] with accuracy exceeding 90%.

In this chapter, we would like to show how the 2D blastula model can exhibit various discrete nonuniform patterns and how these patterns can be quenched by varying kinetic parameters, transport parameters, and the inflow of the substrate.

Advertisement

2. Methods

2.1 Theoretical setup

To model a small multicellular organism of unspecialized cells having only one layer of cells, a blastula, we chose a ring of N = 20 equivalent cells representing the cross section of a blastula. This theoretical setup, see Figure 1, follows the morphogenesis work of A. Turing [16]. The only difference in our case is that both activator and inhibitor are transported at the same rate, as proposed in work by Vastano et al. [23, 24]. Under such circumstances, both nonuniform patterns and uniform oscillations exist simultaneously. If initially, all cells are in the same state, uniform, fully synchronized oscillations appear. In order to obtain a nonuniform stationary pattern, the initial state should be nonhomogeneous. Alternatively, a nonuniform stationary pattern can be obtained by a nonhomogeneous perturbation of synchronized oscillations. This is schematically shown in Figure 1, where in Figure 1a) the 2D blastula model operates in a regime of uniform oscillations, and it acts as 20 synchronized cells. After it is perturbed, it goes into a regime of a discrete nonuniform pattern (discrete Turing pattern) (see Figure 1b). Afterward, it can be switched into other discrete nonuniform patterns (see Results) by an additional perturbation, or it can be switched back to uniform oscillations, as shown in Figure 1c.

Figure 1.

The array of N = 20 coupled cells taken as a 2D blastula model: a) regime of uniform oscillations, b) regime of discrete nonuniform patterns, and c) regime of uniform oscillations again. Black arrows represent the reagent exchange between cells. Each cell is numbered. Blue arrows represent perturbation leading to a transition to discrete nonuniform pattern from uniform oscillations and vice versa. The applied perturbation is shown in Figure 3.

2.2 Model of glycolysis

Our 2D blastula model consists of 20 coupled cells in which chemical reactions occur. A separated cell can be modeled using the kinetic equation:

dZdt=FZP,E1

where vector Z = (x1,…,xm) represents concentrations of reagents involved and vector P = (p1,…,pl) represents the values of l parameters that influence the reaction rates. The vector function F(Z,P) = (F1(Z,P),…,Fm(Z,P)) represents the reaction rates. For N coupled cells in a ring geometry communicating via diffusive transport, the system of evolution equations is

dZdti=QiZiP=FZiP+diagKdZi+12Zi+Zi1,i=1,,N,E2

and the boundary condition is Zi = Zi + N. Vector Kd = (k1,…,km) represents a vector of the transport rate coefficient for each species.

In the case considered in our study, all species have the same transport rate coefficient in between all cells, and thus kd = k1 = k2 = … = km. To describe reaction kinetics, we apply the core model of glycolysis [42] on Eq. (2), yielding to equations:

FxiyiP=νit+σinhyinMn+yinσMxi1+xi1+yi2L+1+xi21+yi2
FyiyiP=ϕσMxi1+xi1+yi2L+1+xi21+yi2kSyiϕσinhyinMn+yinE3

i = 1,..N,

where xi and yi represent ATP and ADP concentration in the i-th cell, respectively. The function νi(t) describes the ATP inflow rate, ϕ if the ratio of dissociation constants of ATP to ADP, n is Hill coefficient, L is allosteric constant, ks is coefficient of degradation of ADP, M is Michaelis constant, σM is autocatalytic rate coefficient, and finally σinh is inhibition rate coefficient. In simulations, we used the values of parameters proposed by Moran and Goldbeter [42] n = 4, L = 106, ks = 0.06 s−1, ϕ = 1. The values of other parameters determining reaction kinetics, σinh, kd, and νi(t), were modified in our simulations. The value of σM = 100 s−1 was used in all simulations.

2.3 Varying parameters to set up the 2D blastula model

Our 2D blastula model has to have carefully chosen parameters, so it can exhibit the coexistence of discrete patterns and nonuniform oscillations. The core model of glycolysis for one cell shows uniform oscillations, stable uniform stationary state, birhythmicity, and hard excitation [42] under fixed autocatalytic rate coefficient σM = 10s−1, varied inhibition rate coefficient σinh ∈ [0 s−1, 4 s−1], and varied ATP inflow rate ν ∈ [0 s−1, 2 s−1]. Since we created a model of arrays of coupled cells with the core model of glycolysis, we had to add transport in between cells originally described by kATP as the transport rate coefficient for ATP and kADP as the transport rate coefficient of ADP. Further studies showed if we increase σM > 80s−1, we can have an equal transport rate coefficient kd for both ATP and ADP. The increment of σM and σinh can be realized experimentally by increasing temperature [43], increasing pH, by the addition of hydrocarbonates [44], or addition of other metabolites [45]. The value of kd ∈ [300;6000] s−1 for whole cellular surface [46]. For intercellular communication, only part of the surface is used.

2.4 Solution diagram

For the analysis of stability and location of bifurcations of stationary states and for simulations of the system, we used the program CONT [47, 48]. A solution diagram was obtained as the result of a one-parameter continuation. The solution diagram for ADP concentration in the first cell as a function of σinh and fixed σM = 100 s−1, kd = 0.1 s−1 is shown in Figure 2. The system has a stable uniform stationary state from σinh ∈ [0 s−1, 17.8 s−1] and σinh ∈ [95.3 s−1, 100 s−1] marked by a solid red line. In between the stable uniform stationary state, there is a region of stable uniform oscillations that occurred via subcritical Hopf bifurcation from the left side and supercritical Hopf bifurcation from the right side, shown by black curves of minima and maxima of concentration of ADP. Inside this region of stable uniform oscillations, there are multiple branches of discrete Turing patterns that occurred from branch point bifurcations marked by a blue dashed line. These patterns are secondarily stabilized by supercritical Hopf bifurcations and, therefore, they cannot occur spontaneously just by varying σinh or kd. The stable discrete nonuniform patterns are marked by solid red curves. The solution diagram only shows possibilities of stationary concentrations of ADP in the first cell. All discrete nonuniform patterns shown simultaneously in all 20 cells are discussed in the Results section.

Figure 2.

Solution diagram of ADP in the first cell for σM = 100 s−1, kd = 0.1 s−1, ν = 1.84 s−1. Red solid curve – Stable stationary state, blue dashed curve – Unstable stationary state, solid black curve – Stable minima and maxima of glycolytic oscillations, dashed black curve – Unstable minima and maxima of glycolytic oscillations, Hopf bifurcation points, and branch points are omitted for the purpose of readability of the solution diagram.

Advertisement

3. Results

Our goal is to provide a method for assessing organism development toxicity. To create the method, we need to map all possible patterns inside the 2D blastula model. For the purpose of machine reading and loading the resulting discrete patterns, we are assigning a letter to certain ranges of stationary concentrations of ADP, Table 1.

LetterRange of stationary concentration of ADP
A[1, 12]
B(12, 28)
C[28, 36]
D(36, 47)
E[47, 100]

Table 1.

Assignation of a letter to a stationary concentration of ADP.

A uniform stationary state has thus pattern C20 with its dimensionless stationary concentration of ADP = 30.6667.

The 2D blastula model can operate in a regime of uniform oscillations, uniform stationary state, and according to Figure 2 at least 12 discrete nonuniform patterns are possible for kd = 0.1 s−1. The analysis of the occurrence of the discrete nonuniform patterns was done by varying the input concentration of ATP by trying all the values within range νi(t) ∈ [0 s−1, 4 s−1] for each cell and also for all combinations of 20 cells.

The same analysis can be done for any chemical input we define as a metabolic pathway inside each cell. The guide for our analysis under constant kd is the solution diagram in Figure 2, or it could be a bifurcation diagram in the whole parameter plane of kd and σinh. A simulation of a transition between uniform oscillations and discrete nonuniform patterns on six randomly chosen cells is shown in Figure 3. The figure is divided into four subfigures illustrating the perturbations and concentrations in the time interval [0 s, 4000 s]. Figure 3a) shows the values of ν1(t), ν6(t), and ν11(t). Figure 3b) shows the values of ν2(t), ν10(t), and ν20(t). Figure 3c) shows the time-dependent concentration of ADP at cell 1, cell 6, and cell 11. Figure 3d) shows the concentration of ADP in cell 2, cell 10, and cell 20. The initial concentrations of both ATP and ADP in every cell are set at value 10, which leads to uniform oscillations. The system oscillates until the first perturbation is applied at 500 seconds. We have decided to use only two types of perturbation values, and that is νi(t) = 3.84 s−1 and no flow of substrate, νi(t) = 0 s−1, while outside perturbation times, the flow of substrate is set according to our previous studies as νi(t) = 1.84 s−1 [25]. The first perturbation is νi(tk= 3.84 s−1 for i = {2,3,6,7,10,11,14,15,18,19} and νj(tk= 0 s−1 for j = {1,4,5,8,9,12,13,16,17,20} applied in the time interval [500 s, 600 s]. It leads to a transition from uniform oscillation to nonuniform pattern type AE2A2E2A2E2A2E2A2E2A. The second perturbation: νi(tl= 0 s−1 for i = {2,3,6,7,10,11,14,15,18,19} and νj(tl= 3.84 s−1 for j = {1,4,5,8,9,12,13,16,17,20} is applied in the time interval [1500s, 1600s]. This perturbation leads to a transition between two nonuniform patterns from type AE2A2E2A2E2A2E2A2E2A to a type EA2E2A2E2A2E2A2E2A2E. The third perturbation has the same νi(t) values as the first perturbation, and it is applied in the time interval [1900s, 2000s]. It leads to a transition between two nonuniform stationary patterns from the type EA2E2A2E2A2E2A2E2A2E to the type AE2A2E2A2E2A2E2A2E2A. All the patterns in this example are the same type, just rotated by two cells toward each other. The fourth perturbation is defined as νi(tz= 3.84 s−1 for i = {1,2,3,6,7,10,11,14,15,18,19} and νj(tz= 0 s−1 for j = {4,5,8,9,12,13,16,17,20}, and it is applied in the time interval [3000 s, 3100 s]. The symmetry of νi(t) values are selected such that they do not correspond to any stationary pattern for this parameter region. The application of such perturbation leads to a transition from pattern type AE2A2E2A2E2A2E2A2E2A back to uniform oscillation. The uniform oscillations prevail until the end of the simulation at 4000 seconds.

Figure 3.

Dynamic simulation of transitions in between uniform oscillations and nonuniform discrete patterns, σM = 100 s−1, σinh = 35 s−1, kd = 0.1 s−1: a) red curve represents ν1, the green curve represents ν11, the blue dashed curve represents ν6; b) red curve represents ν2, the green curve represents ν20, and the blue dashed curve represents ν10; c) ADP concentrations in selected cells: Red curve – Cell #1, green curve – Cell #11, blue curve – Cell #6 d) ADP concentrations in selected cells: Red curve – Cell #2, green curve – Cell #20, blue curve – Cell #10.

To assess all accessible discrete nonuniform patterns, we used initial conditions taken from the solution diagram in Figure 2, or in case any pattern occurs for different parameters than those used in Figure 2, we took the initial conditions from the whole parameter space we have analyzed.

The first set of discrete nonuniform patterns is presented in Figure 4. The patterns in Figure 4a and b) show concentration profile E2A3E2A3E2A3E2A3. They occur for σinh = 60 s−1 and kd = 0.1 s−1. We can also describe the pattern using wavenumber, defined as the number of wavelengths per the distance of 20 cells. For the case of Figure 4a and b), the wavenumber is 4 for pattern E2A3. The patterns in Figure 4c and d) are different as they represent C4DC4DC4DC4D with wavenumber 4. They exist for σinh = 26 s−1, kd = 0.1 s−1. The nonuniform pattern in Figure 4e) has pattern A2BCD4CBA2BCD4CB, and it occurs for σinh = 35 s−1, kd = 0.55 s−1. The wavenumber of it is 2. This pattern is present in five rotations positions within the 2D blastula model. The remaining five rotation positions of this pattern are found for σinh = 60 s−1 and kd = 0.55 s−1, but we decided not to show it as the pattern A2BCD4CBA2BCD4CB is the same.

Figure 4.

Discrete patterns of the stationary concentrations of ADP in 20 cells of the artificial blastula arranged in the time interval [500 s, 1000 s] for various sets of parameters σM = 100 s−1. a) σinh = 60 s−1, kd = 0.1 s−1, b) σinh = 60 s−1, kd = 0.1 s−1, c) σinh = 26 s−1, kd = 0.1 s−1, d) σinh = 26 s−1, kd = 0.1 s−1, and e) σinh = 35 s−1, kd = 0.55 s−1. Each cell is numbered by a green number. The blue dashed line corresponds to ADP concentration in an unstable uniform stationary state. The red lines represent stationary concentrations of ADP in each cell.

The second set of discrete nonuniform patterns is specific. Figure 5 shows nine discrete patterns which coexist simultaneously for σinh = 35 s−1 and kd = 0.1 s−1 and are not all just rotations of the same pattern. The patterns in Figure 5a–d) show concentration profile A2E2A2E2A2E2A2E2A2E2, which has wavenumber 5. The patterns in Figure 5e–g) appear to be somehow connected by their concentration profile values, but after careful examination, they are three different patterns. The pattern in Figure 5e) has the concentration code D2CAD2ACD2A2E2A2E2A2, the pattern in Figure 5f) has the concentration code DEA2E2A2EDACDEA2EDCA, and the Figure 5g) has the concentration code D3A2E2A2D3ADEA2EDA. The wavenumber of all three patterns is 1. These three patterns are highly irregular compared to other patterns, and they may be responsible for either unwanted development mutation or are necessary to create some type of polarity in the organism for development beyond blastula and gastrula. The pattern in Figure 5h and i) has concentration profile A2D3A2D3A2D3A2D3 and the wavenumber 4.

Figure 5.

Discrete patterns of the stationary concentrations of ADP in 20 cells of the artificial blastula arranged in the time interval [500 s, 1000 s] for σM = 100 s−1, σinh = 35 s−1, kd = 0.1 s−1 each cell is numbered by a green number. Notation is shown in Figure 4.

Advertisement

4. Discussion

The coexistence of multiple patterns for the same set of parameters illustrated in Figure 5 allows for multiple morphogenetic results. Our analysis of the model has shown there is not a simple pattern with one minimum of concentration of ADP and one maximum of concentration of ADP, and with the wavenumber 1, which we would expect to start gastrulation. There are, however, three patterns with wavenumber 1, where we can observe certain concentration profiles with both maxima and minima.

This might be the pattern, which could start gastrulation as a result of a concentration profile in 10 cells altogether or all 20 cells. In these cases, we expect the gastrulation to start at D2 in between …CAD2AC… or it can start at A2 in between …DEA2ED… or finally it can occur at A2 in between …DEA2ED… We can notice the last two cases have the same structure. This might be either a mechanism to ensure that two patterns can lead to gastrulation or our universal 2D blastula model incorporating both gastrulation and limb development to multiple organisms with a tail and without a tail. Such organism would have only one axis of body symmetry. The concentration profiles show that a variety of perturbations can lead to the same pattern type, just rotated by a few cells. Since our model can also serve as a model for developing organism or may set prepatterns for future limb development creating multiple axes of body symmetry. The importance of growing legs on a position shifted by π/10, while all limbs have a constant position to each other, is insignificant. What could be a problem if the perturbation switches the pattern from wavenumber 4 to wavenumber 5 or vice versa. A chemical causing this is toxic toward morphogenetic development, because the organism will grow a tail when it is not supposed to or vice versa. If applied to a growing palm, a perturbing chemical will influence a number of digits.

For the purpose of a method of obtaining chemoinformatics toxicity information in living developing organisms, our 2D blastula model serves only as a skeleton hybrid/artificial organism since it only works with ATP and ADP and takes into account only the anaerobic part of glycolysis via kinetic parameters. For the purpose of modeling the pattern behavior of cancer cells [49] due to the Warburg effect [50], it can be extended toward the full model of anaerobic glycolysis using the model proposed by Hynne et al. [51, 52]. If we want to assess the pattern toxicity toward cells with a whole glycolytic reaction chain, including aerobic part of glycolysis [7], it is possible to incorporate artificial mitochondria [53, 54].

Our long-term goal is the creation of an array [55] of artificial cells [56], which could simulate behavior of blastula or even start a shape development, depending on the capabilities of the artificial cell.

Advertisement

5. Conclusions

We have performed an analysis of the stability and bifurcation of stationary states for the 2D model of blastula consisting of 20 coupled cells. The chemistry of each cell is described by the two-variable model of glycolysis (cf. Eq. (3)). The system shows a remarkable number of discrete Turing patterns and can be used to model different phenomena.

One of them is the application of coupled cells as a chemical memory unit. Different patterns can be used for different symbol coding. Figure 2 suggests a straightforward method of switching between discrete patterns. In order to optimize memory, an extensive study on the best strategy for switching between patterns is necessary.

The introduced model can also serve as an artificial living organism for studies on developmental toxicity. We have chosen a parameter plane of σinh or kd while having constant σM = 100 s−1. Using dynamic simulation in program CONT by varying input function νi(t) ∈ [0 s−1, 4 s−1] for each and every i-th cell for time 100 seconds, we have found and classified eight discrete patterns. Each pattern has been classified based on its concentration value by the letters ABCDE. We have found five patterns, which are just rotated within the 2D blastula and therefore can serve for future development or limbs as a prepattern, specifically E2A3 with wavenumber 4, C4D with wavenumber 4, A2BCD4CB with wavenumber 2, A2E2 with wavenumber 5, and A2D3 with wavenumber 4. These patterns can lead either to the organism growing four limbs or five limbs or growing four or five fingers on the limb. We have also found three patterns, which resent a discrete “wave” [cf. Figure 5e–g]. However, these patterns have mirror symmetry and therefore cannot have wavenumber 2. It opens an opportunity for gastrulation, or it can lead to a prepattern for later body symmetry development with one axis. These patterns are specifically D2CAD2ACD2A2E2A2E2A2; DEA2E2A2EDACDEA2EDCA, and D3A2E2A2D3ADEA2EDA. In these cases, we expect the gastrulation to start at D2 in between …CAD2AC… or it can start at A2 in between …DEA2ED… or finally it can occur at A2 in between …DEA2ED…

Our analysis of pattern occurrence inside artificial 2D blastula should generalize the current cheminformatics methods to include toxicological databases toward the potential development of living organisms while not actually harming animals. It can also extend cheminformatics databases about toxicity toward cancer cells, since the model incorporates anaerobic glycolysis and, therefore, can describe the Warburg effect. The model currently works with ATP and ADP but can be extended toward all species taking part in glycolysis. It can also incorporate artificial mitochondria to work with aerobic glycolysis.

Advertisement

Acknowledgments

This publication is part of a project that has received funding from the European Union’s Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement No. 847413.

Scientific work published as part of an international co-financed project founded from the program of the Minister of Science and Higher Education entitled “PMW” in the years 2020–2024; agreement no. 5005/H2020-MSCA-COFUND/2019/2.

Advertisement

Conflict of interest

“The authors declare no conflict of interest.”

Advertisement

Data availability

All numerical data and the model are stored in the repository at: https://doi.org/10.18150/JKFBUW

References

  1. 1. Tongman S, Chanama S, Chanama M, Plaimas K, Lursinsap C. Metabolic pathway synthesis based on predicting compound transformable pairs by using neural classifiers with imbalanced data handling. Expert Systems with Applications. 2017;88:45-57. DOI: 10.1016/j.eswa.2017.06.026
  2. 2. Medina-Franco JL, Thomas Caulfield T. Advances in the computational development of DNA methyltransferase inhibitors. Drug Discovery Today. 2011;16:418-425. DOI: 10.1016/j.drudis.2011.02.003
  3. 3. Enoch SJ, Hasarova Z, Cronin MTD, Bridgwood K, Rao S, Kluxen FM, et al. Sub-structure-based category formation for the prioritisation of genotoxicity hazard assessment for pesticide residues: Sulphonyl ureas. Regulatory Toxicology and Pharmacology. 2022;129:105115. DOI: 10.1016/j.yrtph.2022.105115
  4. 4. Rizvi NF, Santa Maria JP, Nahvi A, Klappenbach J, Klein DJ, Curran PJ, et al. Targeting RNA with small molecules: Identification of selective, RNA-binding small molecules occupying drug-like chemical space. SLAS Discovery. 2020;25:384-396. DOI: 10.1177/2472555219885373
  5. 5. Manigrasso J, Marcia M, De Vivo M. Computer-aided design of RNA-targeted small molecules: A growing need in drug discovery. Chem. 2021;7:2965-2988. DOI: 10.1016/j.chempr.2021.05.021
  6. 6. Ansari MHR, Saher S, Parveen R, Khan W, Khan IA, Ahmad S. Role of gut microbiota metabolism and biotransformation on dietary natural products to human health implications with special reference to biochemoinformatics approach. Journal of Traditional and Complementary Medicine. 2022. DOI: 10.1016/j.jtcme.2022.03.005 [In press]
  7. 7. Wills LP. The use of high-throughput screening techniques to evaluate mitochondrial toxicity. Toxicology. 2017;391:34-41. DOI: 10.1016/j.tox.2017.07.020
  8. 8. Lo YC, Senese S, France B, Gholkar AA, Damoiseaux R, Torres JZ. Computational cell cycle profiling of Cancer cells for prioritizing FDA-approved drugs with repurposing potential. Scientific Reports. 2017;7:2045-2322. DOI: 10.1038/s41598-017-11508-2
  9. 9. Amo del EM, Rimpelä AK, Heikkinen E, Kari OK, Ramsay E, Lajunen T, et al. Pharmacokinetic aspects of retinal drug delivery. Progress in Retinal and Eye Research. 2017;57:134-185. DOI: 0.1016/j.preteyeres.2016.12.001
  10. 10. Flores-Carrillo P, Velázquez-López JM, Aguayo-Ortiz R, Hernández-Campos A, Trejo-Soto PJ, Yépez-Mulia L, et al. Synthesis, antiprotozoal activity, and chemoinformatic analysis of 2-(methylthio)-1H-benzimidazole-5-carboxamide derivatives: Identification of new selective giardicidal and trichomonicidal compounds.vEuropean. Journal of Medicinal Chemistry. 2017;137:211-220. DOI: 10.1016/j.ejmech.2017.05.058
  11. 11. Ghiano DG, Recio-Balsells A, Bortolotti A, Defelipe LA, Turjanski A, Morbidoni HR, et al. New one-pot synthesis of anti-tuberculosis compounds inspired on isoniazid. European Journal of Medicinal Chemistry. 2020;208:112699. DOI: 10.1016/j.ejmech.2020.112699
  12. 12. Gajewicz-Skretna A, Gromelski M, Wyrzykowska E, Furuhama A, Yamamoto H, Suzuki N. Aquatic toxicity (Pre)screening strategy for structurally diverse chemicals: Global or local classification tree models? Ecotoxicology and Environmental Safety. 2021;208:111738. DOI: 10.1016/j.ecoenv.2020.111738
  13. 13. Chan HCS, Shan H, Dahoun T, Vogel H, Yuan S. Advancing drug discovery via artificial intelligence. Trends in Pharmacological Sciences. 2019;40:592-604. DOI: 10.1016/j.tips.2019.06.004
  14. 14. Sansare S, Duran T, Mohammadiarani H, Goyal M, Yenduri G, Costa A, et al. Artificial neural networks in tandem with molecular descriptors as predictive tools for continuous liposome manufacturing. International Journal of Pharmaceutics. 2021;603:120713. DOI: 10.1016/j.ijpharm.2021.120713
  15. 15. Cai D, van Rijsbergen CJ. Learning semantic relatedness from term discrimination information. Expert Systems with Applications. 2009;36:1860-1875. DOI: 10.1016/j.eswa.2007.12.072
  16. 16. Turing A. The chemical basis of morphogenesis. Philosophical Transactions on Royal Society London B. 1952;237:37-72. DOI: 10.1098/rstb.1952.0012
  17. 17. Kozubowski L, Saito K, Johnson JM, Howell AS, Zyla TR, Lew DJ. Symmetry-breaking polarization driven by a Cdc42p GEF-PAK complex. Current Biology. 2008;18:1719-1726. DOI: 10.1016/j.cub.2008.09.060
  18. 18. Castro-e-Silva A, Bernardes AT. Gastrulation as a self-organized symmetry breaking process. Physica A: Statistical Mechanics and its Applications. 2005;352:535-546. DOI: 10.1016/j.physa.2004.10.039
  19. 19. Murray JD. Mathematical Biology II: Spatial Models and Biomedical Applications. Third ed. USA: Springer Science+Business, Media, LLC; 2003. DOI: 10.1007/b98869
  20. 20. Bagudu A, Kraemer C, Germann P, Menshykau D, Iber D. Digit patterning during limb development as a result of the BMP-receptor interaction. Scientific Reports. 2012;991:1-13. DOI: 10.1038/srep00991
  21. 21. Meinhardt M, Gierer A. Application of a theory of biological pattern formation based on lateral inhibition. Journal of Cell Science. 1974;15:321-346
  22. 22. Meinhardt M, Gierer A. Pattern formation by local self-activation and lateral inhibition. BioEssays. 2000;22:753-760. DOI: 10.1002/1521-1878(200008)22:8<753::AID-BIES9>3.0.CO;2-Z
  23. 23. Vastano JA, Pearson JE, Horsthemke W, Swinney HL. Chemical pattern formation with equal diffusion coefficients. Physics Letters A. 1987;124:320-324. DOI: 10.1016/0375-9601(87)90019-3
  24. 24. Vastano JA, Pearson JE, Horsthemke W, Swinney LH. Turing patterns in an open reactor. The Journal of Chemical Physics. 1988;88:6175-6181. DOI: 10.1063/1.454456
  25. 25. Muzika F, Schreiber I. Control of Turing patterns and their usage as sensors, memory arrays, and logic gates. The Journal of Chemical Physics. 2013;139:164108. DOI: 10.1063/1.4825379
  26. 26. Muzika F, Schreiberová L, Schreiber I. Discrete Turing patterns in coupled reaction cells in a cyclic array. Reaction Kinetic Mechanism Catalysis. 2016;118:99-114. DOI: 10.1007/s11144-016-1004-y
  27. 27. Asakura K, Konishi R, Nakatani T, Nakano T, Kamata M. Turing pattern formation by the CIMA reaction in a chemical system consisting of quaternary alkyl ammonium cationic groups. The Journal of Physical Chemistry. B. 2011;115:3959-3963. DOI: 10.1021/jp111584u
  28. 28. Horváth J, Szalai I, De Kepper P. Designing stationary reaction–diffusion patterns in pH self-activated systems. Accounts of Chemical Research. 2018;51(12):3183-3190. DOI: 10.1021/acs.accounts.8b00441
  29. 29. Castets V, Dulos E, Boissonade J, Kepper PD. Experimental evidence of a sustained standing Turing-type nonequilibrium chemical pattern. Physical Review Letters. 1990;64:2953-2956. DOI: 10.1103/PhysRevLett.64.2953
  30. 30. Rudovics B, Barillot E, Davies PW, Dulos E, Boissonade J, Kepper PD. Experimental studies and quantitative Modeling of Turing patterns in the (chlorine dioxide, iodine, malonic acid) reaction. The Journal of Physical Chemistry. A. 1999;103:1790-1800. DOI: 10.1021/jp983210v
  31. 31. Sanz-Anchelergues A, Zhabotinsky AM, Epstein IR, Muñuzuri AP. Turing pattern formation induced by spatially correlated noise. Physical Review E. 2001;63:056124. DOI: 10.1103/PhysRevE.63.056124
  32. 32. Berenstein IB, Dolník M, Yang L, Zhabotinsky AM, Epstein IR. Turing pattern formation in a two-layer system: Superposition and superlattice patterns. Physical Review E. 2004;70:046219. DOI: 10.1103/PhysRevE.70.046219
  33. 33. Bar-Eli K. Coupling of chemical oscillators. The Journal of Physical Chemistry. 1984;88:3616-3622. DOI: 10.1021/j150660a048
  34. 34. Bar-Eli K, Reuveni S. Stable stationary states of coupled chemical oscillators. Experimental evidence. The Journal of Physical Chemistry. 1985;89:1329-1330. DOI: 10.1021/j100254a002
  35. 35. Dolník M, Berenstein I, Zhabotinsky AM, Epstein IR. Spatial periodic forcing of Turing structures. Physical Review Letters. 2001;87:238301. DOI: 10.1103/PhysRevLett.87.238301
  36. 36. Dolník M, Marek M. Extinction of oscillations in forced and coupled reaction cells. The Journal of Physical Chemistry. 1988;92:2452-2455. DOI: 10.1021/j100320a014
  37. 37. Crowley MF, Epstein IR. Experimental and theoretical studies of a coupled chemical oscillator: Phase death, multistability, and in-phase and out-of-phase entrainment. The Journal of Physical Chemistry. 1989;93:2496-2502. DOI: 10.1021/j100343a052
  38. 38. Yoshimoto M, Yoshikawa K, Mori Y. Coupling among three chemical oscillators: Synchronization, phase death, and frustration. Physical Review E. 1993;47:864-874. DOI: 10.1103/PhysRevE.47.864
  39. 39. Muzika F, Valent I. Irregular oscillations and symmetry breaking in a ring of coupled cells with yeast extract. 2022; DOI: 10.18150/GB24WE
  40. 40. Górecki J, Bose A. Computing with networks of chemical oscillators and its application for schizophrenia diagnosis. Frontier in Chemical Science. 2022;10:848685. DOI: 10.3389/fchem.2022.848685
  41. 41. Górecki J, Bose A. How does a simple network of chemical oscillators see the Japanese flag? Frontiers in Chemistry. 2020;09:580703. DOI: 10.3389/fchem.2020.580703
  42. 42. Goldbeter A, Moran F. Onset of birhytmicity in a regulated biochemical system. Biophysical Chemistry. 1984;20:149-156. DOI: 10.1016/0301-4622(84)80014-9
  43. 43. Mair T, Warnke C, Tsuji K, Müller SC. Control of glycolytic oscillations by temperature. Biophysical Journal. 2005;88:639-646. DOI: 10.1529/biophysj.104.043398
  44. 44. Hereng TH, Elgstøen KBP, Eide L, Rosendal KR, Skålhegg BS. Serum albumin and HCO3- regulate separate pools of ATP in human spermatozoa. Human Reproduction. 2014;29:918-930. DOI: 10.1093/humrep/deu028
  45. 45. Mediavilla D, Metón I, Baanate IV. Purification and kinetic properties of 6-phosphofructo-1-kinase from gilthead sea bream muscle. Biochimica et Biophysica Acta. 2007;1770:706-715. DOI: 10.1016/j.bbagen.2006.11.014
  46. 46. Nazarea AD. Spatiotemporal pattern formation in thin layers and membranes: Critical focal size. Proceedings of the National Academy Science. 1974;71:3751-3753. DOI: 10.1073/pnas.75.9.4313
  47. 47. Kubíček M, Marek M. Computational Methods In Bifurcation Theory and Dissipative Structures. New York, NY: Springer Verlag; 1983. p. 243. DOI: 10.1007/978-3-642-85957-1
  48. 48. Kohout M, Schreiber I, Marek M. A computational tool for nonlinear dynamical and bifurcation analysis of chemical engineering problems. Computer Chemical Engineering. 2002;26:517-527. DOI: 10.1016/S0098-1354(01)00783-9
  49. 49. Feng L, Reynisdóttir I, Reynisson J. The effect of PLC-γ2 inhibitors on the growth of human tumour cells. European Journal of Medicinal Chemistry. 2012;54:463-469. DOI: 10.1016/j.ejmech.2012.05.029
  50. 50. Warburg O. On the origin of cancer cells. Science. 1956;123(3191):309-314. DOI: 10.1126/science.123.3191.309
  51. 51. Hynne F, Danø S, Sørensen PG. Full-scale model of glycolysis in Saccharomyces cerevisiae. Biophysical Chemistry. 2001;94:121-163. DOI: 10.1016/S0301-4622(01)00229-0
  52. 52. Herdeen van JH, Wortel MT, Bruggeman FJ, Heijnen JJ, Bollen YJM, Planqué R, et al. Lost in transition: Start-up of glycolysis yields subpopulations of nongrowing cells. Science. 2014;343:1245114. DOI: 10.1126/science.1245114
  53. 53. Kumar S, Karmacharya M, Michael IJ, Choi Y, Kim J, Kim IU, et al. Programmed exosome fusion for energy generation in living cells. Nature Catalysis. 2021;4:763-774. DOI: 1038/s41929-021-00669-z
  54. 54. Thakur A, Johnson A, Jacobs E, Zhang K, Chen J, Wei Z, et al. Energy sources for exosome communication in a Cancer microenvironment. Cancers. 2022;14(7):1698. DOI: 10.3390/cancers14071698
  55. 55. Li DY, Zhou ZH, Yu YL, Deng NN. Microfluidic construction of cytoskeleton-like hydrogel matrix for stabilizing artificial cells. Chemical Engineering Science. 2022;264:118186. DOI: 10.1016/j.ces.2022.118186
  56. 56. Herianto S, Chien PJ, Ho JA, Tu HL. Liposome-based artificial cells: From gene expression to reconstitution of cellular functions and phenotypes. Biomaterials Advances. 2022;142:213156. DOI: 10.1016/j.bioadv.2022.213156

Written By

František Muzika and Jerzy Górecki

Submitted: 20 January 2023 Reviewed: 23 January 2023 Published: 04 March 2023