Cellular “bauplans”: Evolving unicellular forms by means of Julia sets and Pickover biomorphs

BioSystems 98 (2009) 19–30

Contents lists available at ScienceDirect

BioSystems journal homepage: www.elsevier.com/locate/biosystems

Cellular “bauplans”: Evolving unicellular forms by means of Julia sets and Pickover biomorphs Nelly Selem Mojica a , Jorge Navarro b , Pedro C. Marijuán b , Rafael Lahoz-Beltra a,∗ a b

Department of Applied Mathematics, Faculty of Biological Sciences, Complutense University of Madrid, Madrid 28040, Spain Instituto Aragonés de Ciencias de la Salud, Zaragoza 50009, Spain

a r t i c l e

i n f o

Article history: Received 27 April 2009 Received in revised form 29 June 2009 Accepted 2 July 2009 Keywords: Cellular bauplans Pickover biomorphs Morphogenetic field Julia set Evolving fractal Cytoskeletal mechanical forces Organismic form

a b s t r a c t The universe of cellular forms has received scarce attention by mainstream neo-Darwinian views. The possibility that a fundamental trait of biological order may consist upon, or be guided by, developmental processes not completely amenable to natural selection was more akin to previous epochs of biological thought, i.e. the “bauplan” discussion. Thirty years ago, however, Lynn and Tucker studied the biological mechanisms responsible for defining organelles position inside cells. The fact that differentiated structures performing a specific function within the eukaryotic cell (i.e. mitochondrion, vacuole, or chloroplast) were occupying specific positions in the protoplasm was the observational and experimental support of the ‘morphogenetic field’ notion at the cellular level. In the present paper we study the morphogenetic field evolution yielding from an initial population of undifferentiated cells to diversified unicellular organisms as well as specialized eukaryotic cell types. The cells are represented as Julia sets and Pickover biomorphs, simulating the effect of Darwinian natural selection with a simple genetic algorithm. The morphogenetic field “defines” the locations where cells are differentiated or sub-cellular components (or organelles) become organized. It may be realized by different possibilities, one of them by diffusing chemicals along the Turing model. We found that Pickover cells show a higher diversity of size and form than those populations evolved as Julia sets. Another novelty is the way that cellular organelles and cell nucleus fill in the cell, always in dependence on the previous cell definition as Julia set or Pickover biomorph. Our findings support the existence of specific attractors representing the functional and stable form of a differentiated cell—genuine cellular bauplans. The configuration of the morphogenetic field is “attracted” towards one or another attractor depending on the environmental influences as modeled by a particular fitness function. The model promotes the classical discussions of D’Arcy Thompson and the more recent views of Waddington, Goodwin and others that consider organisms as dynamical systems that evolve through a ‘master plan’ of transformations, amenable to natural selection. Intriguingly, the model also connects with current developments on mechanobiology, highlighting the informational–developmental role that cytoskeletons may play. © 2009 Elsevier Ireland Ltd. All rights reserved.

1. Introduction In biological thought, the concept of a body plan has derived from two distinct intellectual paths (Raff, 1996). On the one side, Owen had put forward in the 1840s the notion of archetypes to characterize the essential design that underlies the diversity of species within a group. And on the other side, a little bit earlier, von Baer’s developmental studies were concerned on how the body plan (bauplan) is deeply immanent in the process of ontogeny for each species. The possibility that this discussion, also connected with the great controversies Cuvier–Geoffroy of that time, could apply not

∗ Corresponding author. Tel.: +34 913945243; fax: +34 913945051. E-mail address: [email protected] (R. Lahoz-Beltra). 0303-2647/$ – see front matter © 2009 Elsevier Ireland Ltd. All rights reserved. doi:10.1016/j.biosystems.2009.07.002

only to the structures of multicellular organisms but also to individual cells, was not considered except by the pioneering approach of D’Arcy Thompson, a few decades later. His discussion on surface tension, viscosity, “rolling curves”, and stability as applying to the determination of cellular form was seminal; later on authors such as Waddington, Goodwin and others were inspired by it, as we are going to argue.

1.1. Cellular Bauplans and Morphogenetic Fields The extension of the idea of a “body plan” to the organization of the cell itself was due to Lynn and Tucker (1976), about thirty years ago, who first considered the possibility of two different mechanisms as responsible for defining organelles’ position inside cells. Structural positioning and chemical signaling were instruct-

20

N.S. Mojica et al. / BioSystems 98 (2009) 19–30

ing organelles to follow a well-ordered and characteristic pattern in many unicellular organisms, such as ciliates. The fact that differentiated structures within a cell, such as a mitochondrion, vacuole, or chloroplast, that perform a specific function, were occupying specific positions inside different classes of cells was the observational and experimental support of the “morphogenetic field” notion in the cellular realm. Further, these authors proposed that the spatial complexity of the organelle arrays is comparable with the multicellular arrangement of different cell types, thus tissues (i.e. muscle, nerve, etc.) and organs (structures made up of several tissues capable of performing some special biological function) in animals. Frankel (1975) argued that some of the mechanisms underlying the spatial specification of pattern are the same both in ciliates (a class of protozoan) and multicellular animals. However, whereas the genes dedicated to specifying shape seem important in higher organisms (for example the Drosophila fly), there appear to be very few of such genes in microorganisms (Harold, 1990). In consequence, at least in microorganisms or unicellular organisms, morphogenesis would be a matter of large-scale spatial organization arising epigenetically, not being encoded directly in the genome. In agreement with Harold (1990), this morphogenetic idea had been considered before by Waddington, who introduced the notion of epigenetic landscape—the fact that a morphogenetic surface which specific topography depends on all the interacting genes or their mutations and not on a particular gene or mutation. In consequence, the form of unicellular organisms could not be a direct expression of a genetic program. Explanations of this kind have led to contemporary evo-devo biologists to reconsider the notion of morphogenetic field (Raff, 1996), assuming under a theoretical point of view that living organisms either unicellular or pluricellular are organized by fields; in point of fact revisiting molecularly the classical notion of bauplan, nowadays also extended within the cellular level after the pioneering work of Lynn and Tucker (1976). Properly in developmental studies, the notion of morphogenetic field had been introduced in the 1920s by Spemann (1921, 1938), Gurwitsch (1922), and Weiss (1923, 1939) which also used the terms of “developmental” or “embryonic” fields. Formally, as we will discuss later on, a morphogenetic field may be defined as a plane A with points representing the locations where cells are differentiated or sub-cellular components (or organelles) organized, thus being the plane A either a tissue sheet or a whole cell, respectively. According with Harold (1990), such field has special spatial and temporal features with some ‘agency’ acting in a coordinated way over the area A. For instance, the agency may correspond to signaling molecules defining the morphogen gradients introduced by Turing (1952), or to endogenous electric fields (Nuccitelli, 1984, 1988), to ionic currents (Limozin et al., 1997), visco-elastic field tensions (Briere and Goodwin, 1990), osmotic fields (O’Shea, 1988) or even to electromagnetic fields (Chovnjuk et al., 2001). No matter which one of the aforementioned agencies is the plausible one, morphogenetic fields always consider some physical rule that, in the last term, is responsible for the genesis and control of such field. Nowadays, two main different points of view are considered, which are not totally incompatible between them (even more, for some authors they are elegantly intertwined: Basler and Struhl, 1994; Marijuán, 1996). One paradigm considers that A locations ‘hold’ positional information that controls the morphogenetic destiny or the self-organization events experimented by cells or their sub-cellular components. Under this view, originally based on the morphogen theory proposed by Turing (1952), the local interactions between points (i.e. the activator and inhibitor morphogens) govern the final emergent pattern. At present, this view is supported by the notion of ‘positional information’ or ‘embryonic organizer’ proposed by Wolpert (1969), Spemann and Mangold (1924), as well as by some genetic mechanisms proposed to explain the forma-

tion of morphogen gradients (Gierer, 1981; Gierer and Meinhardt, 1972; Meinhardt, 2003). Alternatively, the other developmental view considers that cells located at points of A into the morphogenetic field follow a ‘master plan’ resulting the final whole organism as a consequence of some ‘attractors’ in the morphospace. This latter perspective, mainly proposed by Waddington (1940, 1942, 1957), Goodwin (1994) and D’Arcy Thompson and Bonner (1992), considers that organisms are dynamical systems. In a mathematical sense, the organism may be defined as a manifold M called the ‘organism state’, and a smooth evolution function (t) that for any element of t, the time, maps a point of the phase space back into the phase space. The organism evolves through transformations which are mainly explained by the molecular organization and the mechanical properties of the cytoplasm (i.e. via cytoskeleton and networks of filamentous protein polymers, see Hotani et al., 1992; Goodwin and Trainor, 1985). That cytoskeletons may create a ‘mechanical’ environment guiding cellular development and differentiation is receiving strong empirical support. In a very challenging paper, Ainsworth (2008) brings forward two striking experiments of how a mechanical field of physical forces (i.e. governed by cytoskeleton) could be affecting the cells fate. One of these experiments describes how stem cells growing on soft, medium or rigid matrices resembling their natural environments start developing into neurons, muscle and bone cells. Another experiment is related with the mechanics of the cancerous cell: the extracellular matrix or microenvironment surrounding cells with mutated oncogenes becomes stiffer even before those cells form tumours. We will discuss these challenging views, and how mechanical forces, together with chemical morphogens, could organize the developmental and cell differentiation processes of an organism. 1.2. Modeling the Morphogenetic Field The evolution of the morphogenetic field is a relevant topic in Biology since the ways in which morphogenetic fields can change over time may yield insights into the evolutionary possibilities of unicellular and pluricellular organisms. For instance, Perales-Graván and Lahoz-Beltra (2004) simulated the evolution of the morphogenetic field in the zebra skin pattern. In the aforementioned paper we used a genetic algorithm with chromosomes coding the diffusion distance and the morphogenetic constant field values of two different morphogen molecules, the activator and inhibitor, codifying such values the transition rules of a cellular automata population. This cellular automata model of the morphogenetic field of the zebra skin pattern (Young, 1984) was the discrete version of the reaction–diffusion equations proposed by Turing (1952). Surprisingly we found that during evolution spot patterns appear more often than stripe patterns on the simulated skin of zebras, concluding that the stripe pattern of zebras may be a result of other biological features (i.e. genetic interactions, the Kipling hypothesis, etc.) not included in this model—see Perales-Graván and Lahoz-Beltra (2004). Assuming the ‘positional paradigm’, where A locations contain information governing the morphogenetic events, two principal mechanisms have been proposed to convey positional information in the mathematical models (Kavlock and Setzer, 1996). The most common model is the reaction–diffusion approach based on Turing’s morphogen theory, where pairs of diffusing chemicals interact with each other to form stable patterns (i.e. zebra strip patterns). One example of this model is the aforementioned paper (Perales-Graván and Lahoz-Beltra, 2004) where we simulated the evolution of the morphogenetic field in the zebra skin pattern. In the other approach, “genetic addresses” and gradients of chemicals provide the positional information (Basler and Struhl, 1994;

N.S. Mojica et al. / BioSystems 98 (2009) 19–30

Marijuán, 1996). For instance, based on this approach Levin (1994a) supposes that the chemicals of the genes (similar to ‘proteins’) interact intracellularly producing complicate developing patterns, modeling such patterns with fractal and chaos theory. The author also considers that functions in complex variables can be used to simulate genetic interactions, resulting in position-dependent differentiation. Further, this is shown to be equivalent to computing modified Julia sets, and may produce a very rich set of morphologies. Following the assumptions and restrictions introduced by Levin (1994a), in the present paper we study the morphogenetic field evolution yielding from an initial population of cells to different unicellular and microbial organisms (Perry et al., 2002) as well as specialized eukaryotic cells (Lodish et al., 2000). Both types of cells are represented as Julia sets and Pickover biomorphs (Pickover, 1986, 1990) using to simulate Darwinian natural selection a simple genetic algorithm. The rest of this paper has been organized as follows. After this introductory section, Section 2 presents the model and includes a brief description of the main theoretical notions and definitions related with cellular–biological forms, describing the chromosome coding of complex mathematical functions, plus the evolving biomorphs (unicellular organisms as well as specialized eukaryotic cells) obtained by a simple genetic algorithm with a fitness function which evaluates biomorphs based on thirteen selected features related with the form and physical properties of a biomorph. Section 3 presents the computer simulation experiments, and Section 4 makes further comments on the whole results of the simulation experiments. Finally, Section 5 discusses the possible impact of this work and its future directions.

21

Fig. 1. Pickover biomorphs. Different organic-looking shapes obtained as representation of complex functions built with zz , zn , ez , sin(z), etc., terms for complex z variable (from Pickover, 1990).

(|Re(P n (z))| ≤ ε) ∨ (|Im(P n (z))| ≤ ε),

ε∈R

2. Model Let P(z) denote a function in a complex plane C and z0 (z0 ∈ C) a starting value. Given Pn (z) as the function resulting of the composition of P with itself n-times, then the sequence of points P(z), P2 (z), . . ., Pn (z), etc., given by the recursive family zn+1 = Pn (z) with n varying over the set of natural numbers, defines an orbit O(z0 ) of its starting point z0 (Devaney, 1994). For a given function P(z), the behavior of the orbit O(z0 ) defined by the set or sequence of points Pn (z0 ) depends on the selected z0 value. In complex dynamics the set of points the orbit of which is chaotic is called the Julia set. For a given function P(z), the behavior of the orbit O(z0 ) defined by the set or sequence of points Pn (z0 ) is studied using an escapetime algorithm. An escape-time algorithm terminates an iterating formula when either the size of the positive orbit O(z0 ) exceeds a selected bailout real value ε or an iteration limit ı is reached. Once one of the above conditions is reached then the pixel corresponding to the starting point z0 is colored according to the final or last iteration number. For a given Julia set the classical approach obtains the magnitude of Pn (z), thus |Pn (z)| or Euclidean norm of the orbit O(z0 ): |P n (z)| =



Re(P n (z))2 + Im(P n (z))2

where Re(Pn (z)) and Im(Pn (z)) are the real and imaginary parts of z. In such a case, the convergence criterion of the classical escape-time algorithm is given by:



Re(P n (z))2 + Im(P n (z))2 ≤ ε,

ε∈R

In the computer simulation realm, Pickover (1986, 1990) searched for new sets relaxing the convergence criterion, introducing the following criterion:

which has been called in the present paper as the ‘Pickover criterion’. Considering the similarity between living organisms (i.e. unicellular and microbial organisms, see Perry et al., 2002) and the forms obtained by Pickover, such forms were named as ‘biomorphs’ (Fig. 1). Biomorphs can be considered as a particular case of Julia sets. Thus, to obtain a biomorph the escape-time algorithm terminates an iterating formula if either the real or imaginary component is lesser or equal to a bailout value, then a pixel is set to a particular color plotting the biomorph. Note that using such criterion a lesser number of points escaping to infinitum is obtained than with the Euclidean criterion. Based on the above two criteria, in the present paper the escape-time algorithm was defined as follows: 1. Select in the complex plane C a region A. The morphogenetic field is defined as the plane A with points representing the locations where the sub-cellular components or organelles are organized, being the plane A occupied by a cell. 2. Set a maximum number of iterations ı and a bailout value ε (ε ∈ R). 3. Evaluate Pn (z) into a region A until either the size of the orbit O(z0 ) exceeds a selected bailout value ε, thus |Pn (z)| ≥ ε with n ∈ [0,ı], or until the iteration limit is reached. Once the iteration limit ı is reached we compare: a. Given a Julia set, if |Pı (Z)| < ε then the pixel corresponding to the starting point z0 is colored black. Otherwise, the pixel is colored white. b. Given a ‘Pickover set’, if min{|Re(Pı (z))|,|Im(Pı (z))|} < ε then the pixel corresponding to the starting point z0 is colored black. Otherwise, the pixel is colored white. The above algorithm could be written as the following pseudocode program:

22

N.S. Mojica et al. / BioSystems 98 (2009) 19–30

Note how in the Pickover criterion the divergent points also diverge under the Euclidean criterion. Given a Julia set we could write a similar program with the exception of the convergence criterion:

It is important to note how the region A representing the morphogenetic field assumes two X, Y orthogonally diffusing chemicals (Levin, 1994a). Such chemicals are modeled considering X and Y as the real and imaginary components on the complex plane where a biomorph (a cell and its organelles) is ‘drawn’. Under this assumption a plausible schema is: genes → morphogens → cell differentiation. However, in this paper we extended the meaning of X and Y to other possibilities. For instance, if we assume a mechanical field of physical forces, as Ainsworth (2008) proposes, then X and Y could be related to the expression of cell shape genes, somehow ‘mechanical genes’ such as twist and fog (Dawes-Hoang et al., 2005), dacA-rodA-pbpA (De Pedro, 2004), etc. Shape genes would mediate cellular shape and differentiation based on mechanical forces via cytoskeleton, so that: genes → mechanical genes expression → mechanical forces (cytoskeleton) → cell differentiation. Another scheme to ponder is that morphogens may promote the expression of the mechanical genes that rule cell shape and differentiation via cytoskeleton according to the following scheme: genes → morphogens → mechanical genes expression → mechanical forces (cytoskeleton) → cell differentiation. In both schemata, pattern formation is explained by differential gene activation modelled with X, Y (orthogonal vectors) representing physical forces governed by cytoskeleton. Thus, the Re(z) and Im(z) components should be representing physical forces (i.e. through cytoskeleton) on a region A rather than concentrations of gene products as is assumed in Levin’s morphogenetic field.

2.1. Cellular and Organelle Coloring Algorithm In the present paper we model the living matter within a cell assuming that the living matter is composed by cytosol,

organelles and cell nucleus. These cellular elements are represented by pixels with a bounded orbit, being such orbit shown in black color. In order to model living matter we have separated these black pixels in groups of different colors. We begin by selecting n colors and defining a sequence of natural numbers G0 , G1 , . . . , Gn ∈ N such that 0 < G0 < G1 < · · · < Gn = ı, assigning to each interval [Gi−1 , Gi ] a particular and arbitrary color. Following, in the escape-time algorithm we search for the first n such that |Pn (z)| ≥ ε, defining n as ˛, otherwise ˛ = ı. Once the iterations are finished in the Euclidean escape-time algorithm then the Pickover criterion is used: if min{|Re(Pı (z))|,|Im(Pı (z))|} < ε then the pixel is colored depending on the color assigned to the interval [Gi−1 , Gi ], such that Gi−1 ≤ ˛ < Gi . Thus, ˛ is a measure of the escape velocity to infinitum under the Euclidean norm. Finally, we assigned to each color a ‘cellular component’, representing black, dark yellow, red and grey colors the cellular membrane, cytosol, cell nucleus and cellular environment, respectively. The remainder colors were assigned to ‘organelles’ without a distinction between the particular type of organelle and its biological function (i.e. chloroplast, mitochondrion, endoplasmic reticulum, etc.). Note that colored pixels, excepting orange pixels, represent points escaping to infinitum under the Euclidean norm, but not under the Pickover criterion. The orange color was assigned to represent organelles which points do not escape under both criteria. Thus, given a polynomial function P(z) orange color represents the set of points which border is the Julia set. In the present model (see Table I, Lahoz-Beltra et al., 2008) the white color was assigned to those points that diverge very fast under the Euclidean norm but not under Pickover criterion. The grey color was assigned to those

N.S. Mojica et al. / BioSystems 98 (2009) 19–30

points that diverge under both criteria, thus when |Pn (z)| ≥ ε and min{|Re(Pı (z))|,|Im(Pı (z))|} ≥ ε. 2.2. Chromosome Coding of Complex Functions In this model we define a population of cells represented by a set of chromosomes simulated as strings with genes modeling different features of a complex function. In the genetic algorithm terminology, genes were defined as follows. At each gene, and from left to right, the first and second gene positions are integer values representing the complex function terms:

fi (z)(i=1,2) =

⎧ genei = 0, Z√n ⎪ ⎪ n ⎪ ⎪ genei = 1, Z ⎪ ⎪ ⎪ genei = 2, Z z ⎪ ⎨ genei = 3,

exp(z)

genei = 4, cosh(z) ⎪ ⎪ ⎪ ⎪ genei = 5, sinh(z) ⎪ ⎪ ⎪ ⎪ ⎩ genei = 6, sin(z) genei = 7,

cos(z)

whereas the third gene represents the membrane function fm and the fourth gene the complex function f(z) obtained from f1 (z) and f2 (z) terms as follows:

f (z) =

⎧ gene4 = 0, f1 (z) + f2 (z) + c ⎪ ⎪ ⎪ gene4 = 1, Z.f1 (z) + c ⎪ ⎨ gene4 ⎪ gene4 ⎪ ⎪ ⎪ ⎩ gene4 gene4

= 2, = 3, = 4, = 5,

⎧ gene4 = 0, f1 (z) + f2 (z) + c ⎪ ⎪ ⎪ gene4 = 1, Z · f1 (z) + c ⎪ ⎨ gene4 gene4 ⎪ ⎪ ⎪ ⎪ ⎩ gene4 gene4

= 2, = 3, = 4, = 5,

2.3. Evolving Biomorphs with a Simple Genetic Algorithm Genetic algorithms are stochastic optimization procedures based on Darwinian natural selection (Lahoz-Beltra, 2004, 2008). This optimization technique uses genetic operators, thus procedures inspired in the genetic mechanism observed in populations such as crossover or recombination (combination of two solutions) and mutation (random change of a solution). At present, the application of genetic algorithms to fractals is a common approach in applied and theoretical research (Lutton, 1999; Ashlock and Jamieson, 2007; Vences and Rudomin, 1997). In this paper a simple version of a genetic algorithm (Goldberg, 1989) was applied to evolve a population P(t), modeling P(t) a population of ‘cells’ (or biomorphs). Each cell was represented by a complex function. The genetic algorithm was defined as follows: 1. 2. 3. 4.

Choose initial population P(0). Evaluate the fitness of each cell in the population P(t). Select the best-fitness cells to reproduce. Breed a new generation of cells through crossover and mutation and give birth to offspring. 5. Replace worst ranked part of population with offspring. 6. Until . or by the following pseudo-code program:

(0, Im(z)) · f2 (z) + Z + c f1 (z) + c Z · f2 (z) + Z 2 + c Z · f1 (z)2 + Z + c

The membrane function fm is a computational procedure oriented to polish the biomorph membrane which is defined as follows: if the membrane thickness is high (the number of black pixels is greater than one) then remove all those black pixels without contact with colored pixels rendering a thin membrane around cytoplasm. The fifth gene has a binary value, such that if gene5 = 0 then f(z) remains as is without any possibility of transformation. Otherwise, when gene5 = 1 the complex function f(z) is transformed to a new one (evaluating the function in Z2 instead z if the function is trigonometric). Note that transformation is only carried out if gene4 is equal to 2, 4 or 5 values:

f (z) =

23

(0, Im(z)) · f2 (z) + Z 2 + c f1 (z) + c Z 2 · f2 (z) + Z 3 + c Z · f1 (z)2 + Z 2 + c

Finally, the last four genes set the value of the complex constant c. Sixth and seventh genes codify the integer as well as decimal values of the real part of c, whereas the genes number eight and nine codify the integer and decimal values of the imaginary part of the complex constant c. In order to illustrate how a complex function is codified in a chromosome lets go to consider the following example. For instance, let us denote by gene1 = 5, gene2 = 3, gene3 = 1, gene4 = 0, gene5 = 0, gene6 = 4, gene7 = 39, gene8 = 1 and gene9 = −25 a chromosome. Hence, the complex function codified in such chromosome is: f (z) = sinh(z) + exp(z) + (4.39 + 0.75i) Note that in the above expression the complex constant was obtained reading the values for genes 6, 7, 8 and 9 as follows: c = (4 + 0.39) + (1i − 0.25i).

In this model, we define a population of chromosomes simulated as strings with values representing the complex function terms. The genetic algorithm uses two-points recombination and mutation operators with a population size of 200 biomorphs. The algorithm cycles through epochs searching for an optimum biomorph until a maximum of 10 generations is reached, setting the recombination and mutation probabilities to pr = 0.70 and pm = 0.08. Starting with a random population of chromosomes reproduction, recombination and mutation were simulated, obtaining new biomorph generations of equal population size. The initial population of chromosomes was obtained choosing the gene values from a uniform distribution according with the following ranges: gene1 ∈ [0,7], gene2 ∈ [0,7], gene3 ∈ [0,1], gene4 ∈ [0,5], gene5 ∈ [0,1], gene6 ∈ [−9,9], gene7 ∈ [−99,99], gene8 ∈ [−9,9] and gene9 ∈ [−99,99]. 2.3.1. Selection Once the chromosomes are evaluated, we select the mating pool of the next generation using the roulette wheel parents selection algorithm (Davis, 1991; Lahoz-Beltra, 2004). This is a method for implementing reproduction, and thereby Darwinian natural selection, by spinning a roulette wheel that assigns to each chromosome a slot whose arc size is proportional to its fitness value. Of course, other selection schemes are possible such as tournament selection, truncation selection as well as linear and exponential ranking selection (Blickle and Thiele, 1995), however the roulette wheel parents

24

N.S. Mojica et al. / BioSystems 98 (2009) 19–30

selection scheme bears a better resemblance to Darwinian natural selection (Lahoz-Beltra, 2001; Perales-Graván and Lahoz-Beltra, 2004). At each generation, the fitness of each chromosome, thus the degree of achievement of a biomorph was evaluated using the following fairly fitness function: f = prs + pas + pcs + pos + pI + pM + pCm + pCn + pCc + pCo + pCa−p + pCe + pCo−o where prs , pas , pcs , pos , pI , pM , pCm , pCn , pCc , pCo , pCa − p , pCe and pCo − o are evaluation terms which were defined to evaluate in a given biomorph its fitness. Thus, and in the context of genetic algorithms the calculated fitness (see Lahoz-Beltra et al., 2008) represents in each biomorph a measure of the adaptation and the ability to produce offspring in such environment. In consequence, in the present model fitness is a measure of the degree of achievement and survival for a given genotype or complex function. The above evaluation terms were included in the fitness function representing thirteen selected features related with the form and physical properties of biomorphs: the radial symmetry, axial symmetry, cytosol symmetry, organelles symmetry, moment of inertia, total mass, mass of membrane, mass of cell nucleus, mass of cytosol, mass of organelles, membrane interactions and environmental interactions as well as the organelles form complexity, respectively. Next sections define the main features selected and included in the fitness function to evaluate biomorphs. 2.3.1.1. Radial Symmetry. Given a complex plane C a region A has radial symmetry related to a central point p when a rotation T of A results again in the region A, thus T[A] = A. In a biological realm, those biomorphs exhibiting this property will have several cutting planes producing roughly identical pieces, having the biomorph a top and a bottom named as dorsal and ventral surfaces respectively. The algorithm applied to biomorph images was as follows. First, given a biomorph image A then a point p is selected being p the mass center of de image. Secondly, we choose a rotation angle ω = 2k/40 with k in the interval [0,40]. Following, given a ω value we count the number of points (pixels) into the rotated biomorph image or T[A] that coincide with the image or region A before rotation. In the present paper we adopted as practical criterion that is sufficient if T[z] neighborhood contains a point or pixel with the same color of z. Afterwards, we calculate the division between number of points or pixels that coincide and the number of points or pixels that after rotation T[A] remain into the region A. In third place, and assuming that k has been increased one unit at each time, we obtain a list including 40 values such that for a given maximum of the list if the rotational symmetry coefficient Csr is greater than ı (with ı = 0.8) then we enter in the accounts ‘one symmetry’. Finally, once we obtained all the radial symmetries k, thus we found the total number symmetries, a radial symmetry coefficient Cr was calculated as the following average of symmetries:

r

Cr =

C i=1 sr

nmax

where nmax is the total number of maxima. It is important to note that since we are considering cells, thus unicellular organisms and specialized eukaryotic cells, and not multicellular organisms such as plants or animals, then we consider symmetry as the balanced distribution of duplicate organelles instead of the body parts or shapes of animals and plants. 2.3.1.2. Axial Symmetry. Given a complex plane C a region A has axial symmetry when the reflection related to an axis L is invariant. The algorithm applied to biomorph images was similar to the above one but replacing rotations by reflections. First, we select 40 straight lines Lk with slope mk equal to the tangent of ω = 2k/40 (with k

in the interval [0,40]) going the straight lines across of the mass center of a biomorph image or region A. Secondly, we count the total number of axial symmetries and the maximum of axial symmetries (thus, the maximum symmetry coefficient). In the present paper we adopted the following practical criterion. Given an axis L we enter in the accounts as axial symmetry those cases that after reflection the percentage of points or pixels with the same color is greater than ı (with ı = 0.8), being also greater such percentage with regard to the preceding and subsequent reflection axis. 2.3.1.3. Organelles and Cytosol Symmetries. This property is a measure of the unbalance or asymmetry in the distribution of the organelles (or colors) inside a cell or biomorph. Since the cytosol (dark yellow) is the fluid portion of the cytoplasm, outside the organelles, this element has been omitted in such measure. Furthermore, considering that blue color modeling some organelles is the most often (this fact was observed in preliminary simulation experiments) after the dark yellow assigned to cytosol, then the organelles asymmetry has been measured considering only the blue organelles. The algorithm applied to biomorph images was as follows. First, given a biomorph image A then the cytoplasm mass center cmc and the mass center of the blue organelles omc are both calculated, obtaining the two following coefficients |p − cmc | and |p − cmo | which were named as cytosol and organelles asymmetries respectively, being p the mass center of the biomorph. Note that if blue organelles were symmetrically distributed into the cell or biomorph then omc and cmc would be equal between them and equal to the mass center p of the biomorph, thus omc = cmc = p. Otherwise, the above coefficients would be |p − cmc | > 0 and |p − cmo | > 0 such that a greater values higher will be the cytosol asymmetry and organelles symmetry. 2.3.1.4. Inertia. In the present model, inertia has been defined as the tendency of the biomorph ‘body’ to resist rotation, instead of resisting acceleration as is usually defined in Physics. In the biological realm, unicellular organisms like bacteria and protozoan (single-celled eukaryotic organisms, i.e. amoebas, trypanosomes, sporozoans, and paramecia) exhibit cellular appendages specialized for locomotion (i.e. cilium, flagellum). Since inertia is related with mass distribution then biomorph motility and in consequence the food search will be both related with this property. Thus, a biomorph searching for food will change its position and in consequence its mass distribution and inertia being an evolutive advantage for a biomorph to exhibit the radial symmetry. Let mi denote the mass at point i and r the distance from the point i to the mass center p. Given a biomorph image A then the moment of inertia of the biomorph image (IA ) is defined by: IA =



mi r 2

i∈A

2.3.1.5. Mass. This property is defined as biomorph resistance to changes in either the speed or direction of its motion. Assuming that we assign a zero weight to the grey pixels representing the cellular environment (virtually every habitat present in Nature), and a unity weight value to the colored pixels, thus those pixels representing cytosol (dark yellow) and organelles. Given a biomorph image A then the total mass of the biomorph (MA ) has been measured as: MA =

number of colored pixels total number of pixels

being the total number of pixels the sum of colored pixels and grey pixels. Once MA was calculated we proceeded to calculate the mass of the cell nucleus, membrane, cytosol and organelles with regard to the total mass MA of biomorph. At first, biomorphs should contain more cytosol than other cellular structures of the fact that is in such

N.S. Mojica et al. / BioSystems 98 (2009) 19–30

media where metabolism is carried out in biomorphs. Thus, given a biomorph image A then we defined the following mass coefficients. A membrane coefficient Cm representing the mass of membrane respect to the biomorph total mass is given by: Co−c =

Cm =

number of black pixels number of colored pixels − number of black pixels

Cn =

number of red pixels number of colored pixels − number of black pixels

Cc =

number of dark yellow pixels number of colored pixels − number of black pixels

Note that both coefficients are calculated without considering the mass of the biomorph membrane. Finally, the organelles coefficient Co was defined as: number of not red and not dark yellow pixels number of colored pixels − number of black pixels

where the mass of organelles is given by their mass without consider the mass of the biomorph membrane. 2.3.1.6. Membrane Interactions. In the present model different kind of interactions were measured and included in the fitness function. First, an area–perimeter coefficient Ca − p was introduced measuring the amount of real membrane with regard to the non-grey area of the biomorph. It is important to note that a large value of this coefficient indicates that a biomorph is a spore. Thus, a reproductive structure adapted for dispersion and surviving for extended periods of time in unfavorable conditions. For this reason we consider in the model that a ‘spore biomorph’ is not viable at present state until the environment changes to favorable conditions. Spores form part of the life cycles of many plants, algae, fungi and some protozoans. The coefficient was defined as follows: Ca−p =

number of black pixels number of non-grey pixels

Another coefficient, the environmental interactions coefficient Ce , was introduced as a measure of the amount of interactions between the environment and the biomorph membrane, which is given by: Ce =

n x − Imin

n Imax

pixels indicate a high order. However, an intermediate amount of ‘blue pixel–blue pixel’ contacts will be interpreted as a measure of complexity, penalizing those biomorphs with a high or low value of the following coefficient:

number of contacts blue pixel–blue pixel number of contacts blue pixel–blue pixel + number of contacts blue pixel–no blue pixel

whereas for the cell nucleus a coefficient Cn , and for the cytosol a coefficient Cc were both defined representing the mass of cell nucleus and the mass of cytosol with regard to the total mass of biomorph:

Co =

25

n − Imin

where n is the number of pixels of the membrane, x the number x x of neighbors, and Imax and Imin the maximum and minimum number of neighbors considered as possible contacts. Now we assume √ x x = 2n + 2 and Imin = 4 m with m being the smaller square that Imax is bigger than n. Note that high values of this coefficient indicate those biomorphs with a high number of interactions between membrane and the environment. 2.3.1.7. Organelles Complexity. We propose a last coefficient which is a measure of the organelles form complexity. In principle since blue organelles are abundant into biomorphs it is possible to get a measure of order only considering those organelles represented by blue color. Thus a high or low number of contacts between blue

2.3.2. Recombination After a new generation of offspring chromosomes is obtained, a two point crossover occurs, but only between the chromosome segments including the last four genes which codify for the value of the complex constant c. Once two points are randomly selected on the chromosomes, then the segments between the two points are swapped between the chromosomes, rendering two new chromosomes. The decision to perform crossover on a current pair of chromosomes underlies on a Bernoulli trial being the recombination probability pr of 0.70. 2.3.3. Mutation Mutation at a gene was simulated changing at random the value gene, choosing the mutated value from a uniform distribution with similar range that those defined to obtain the initial population of chromosomes. Once again whether or not to change a gene value on a chromosome is decided on the basis of a Bernoulli trial, mutation being a success with a given probability pm (mutation probability) of 0.08.

3. Genetic Algorithm Simulation Experiments Several simulation experiments were carried out assuming in each one that an initial population of cells (or biomorphs), the form of which was plotted from chromosomes, could be evolved in different environments (or under different fitness functions, see Lahoz-Beltra et al., 2008) to one of the following cellular classes: • Class I. This class is defined by ancient cells without cell nucleus and mitochondrion, having a thin membrane. These cells live in a watery environment. • Class II. Is a kind of cells similar to class I excepting that the cells have mitochondrion. • Class III. The cells belonging to this group have small size and a large membrane. • Class IV. This class is defined by small and round cells, with cell nucleus and mitochondrion, having a thin membrane. • Class V. Is a sort of cells similar to class III excepting that the cells have a larger size. • Class VI. The cells belonging to this group are large and round cells, with cell nucleus and organelles. • Class VII. This class is defined by large cells with radial symmetry, having cell nucleus and organelles. • Class VIII. The cells belonging to this group are aligned rectangular cells, with axial symmetry, developing layers or tissues. • Class IX. Is a sort of cells specialized in communication cell-tocell (i.e. neurons). This kind of cell has a membrane with a large surface, this feature promoting a maximum of interactions or ‘synapses’ between a cell and its neighborhood. • Class X. This class is defined by ‘feed cells’ or cells specialized in the absorption of nutrients (i.e. intestinal cells) from environment. Considering their function such cells have axial symmetry and an asymmetrical distribution of organelles into the cell. In con-

26

N.S. Mojica et al. / BioSystems 98 (2009) 19–30

Fig. 2. Evolved Pickover biomorphs. Unicellular ‘organisms’ resembling watery bacteria (class I), aerobic bacteria (class II), bacterial spores (class III), ancient and small protozoans (class IV), protozoan spores (class V), watery protozoans (class VI), protozoans (class VII), and specialized eucariotic animal ‘cells’ bearing a resemblance to epithelial cells (class VIII), neuronal cells (class IX), intestinal cells (class X) and erythrocytes (class XI), were evolved under different simulated ‘environments’ (fitness function) from a random initial population of cells.

sequence, is in the top of cell where is located the food entrance site into the cell. • Class XI. Is a kind of cells specialized in the transport of substances, without cell nucleus (i.e. oxygen is carried by erythrocytes or red blood cells). In consequence, the initial population will evolve to one of the above classes (Fig. 2) depending on the values set up for the evaluation terms of the fitness function (see Lahoz-Beltra et al., 2008). Under a biological realm it is important to note that the proposed classes bear a strong resemblance to the following cellular types. Unicellular or microbial organisms (Perry et al., 2002; Fig. 3a) such as watery bacteria (class I), aerobic bacteria (class II), bacterial

spores (class III), ancient and small protozoans (class IV), protozoan spores (class V), watery protozoans (class VI), protozoans (class VII). Specialized eukaryotic cells (Lodish et al., 2000; Fig. 3b) from multicellular organisms, in particular animals, such as epithelial cells (class VIII), neuronal cells (class IX), intestinal cells (class X) and erythrocytes (class XI), respectively. 4. Simulation Results Computer simulation experiments were carried out using the population size, gene values as well as the recombination and mutation probabilities referred in Section 2.3. In Fig. 4 we show two representative performance graphs of the simulation experiments.

Fig. 3. (a) Prokaryotic cell, i.e. bacterium. (b) Eucariotic cell, showing subcellular components or organelles: (1) nucleolus, (2) nucleus, (3) ribosome, (4) vesicle, (5) rough endoplasmic reticulum, (6) Golgi apparatus, (7) cytoskeleton, (8) smooth endoplasmic reticulum, (9) mitochondria, (10) lysosome, (11) cytoplasm, (12) peroxisome, and (13) centrioles.

N.S. Mojica et al. / BioSystems 98 (2009) 19–30

Fig. 4. Performance graph. (a) Evolutionary dynamics of class IV (ancient and small protozoans) is also displayed by classes III and V. (b) Evolutionary dynamics exhibited by class IX (neuronal cells) is also displayed by classes I, II, VI–VIII, and X.

Performance was measured as the average fitness per generation as is usual in experiments based on genetic algorithms. The evolutionary dynamics of class IV (Fig. 4a) is also displayed by classes III and V. In a similar fashion the dynamics exhibited by class IX (Fig. 4b) is also displayed by classes I, II, VI–VIII, and X. The three first simulation experiments (classes I–III, Fig. 2) were addressed in order to evolve ‘bacteria’ whereas ‘protozoans’ were evolved under class IV to class VII experiments (Fig. 2). In the class I experiment we found that in the last generation the evolved bacterial biomorphs have a thin membrane being cells with a rounded form and small size. Such cells exhibit a symmetry in the colors distribution and those cells including blue organelles evolved the aforementioned organelles very early in the course of evolution, thus during the firsts generations. However, compared these results with the biomorphs evolved under class II experiment, those bacterial cells belonging to class II show a variety of organelles, mainly represented by blue, green and orange colors, compared with almost the absence of organelles in class I. They are cells of small size, thin membrane and radial symmetry. One difference between both classes is that in class I the average fitness value is very close to the maximum fitness, whereas in class II the average fitness does not reaches the maximum value. The results obtained for class III experiments were ‘bacterial spores’. This class of biomorph exhibits a large membrane, being cells ‘more resistant’ to the environment with enough cytosol (dark yellow) inside cell to survive. Spores performance graph shows a sudden fitness jump in the first generation and how the average fitness growths very slowly after the jump. Compared with above bacterial classes, in the experiments evolving protozoan cells at least the 20% of the evolved cells of class

27

IV (Fig. 4a) exhibit a cell nucleus (red) being protozoan biomorphs with a thin membrane and a wide variety of forms as result of the promotion of both, the radial symmetry and axial symmetry. Likewise the bacterial spores of class III, protozoan spores of class V have a large membrane. However the protozoan spores exhibit a higher number of blue organelles, obtaining as a deceitful result that the number of biomorphs with cell nucleus is only about the 15% of the evolved cells. The class VI is composed by evolved cells with a high biodiversity of cell types. Thus, they are cells of big size and colored without any restriction about the type of symmetry. Such protozoans exhibit a thin membrane and, even when natural selection does not promote a particular type of symmetry, there is often a number of cells having axial symmetry. Finally, the protozoan biomorphs belonging to class VII are big cells with cell nucleus and organelles, having same of them radial symmetry. Furthermore, the obtained cells are very similar to the intestinal cells of class X. The last three experiments concern with specialized eukaryotic cells present in animal tissues such as skin, brain, intestine and blood, the cells of which were evolved under classes VIII, IX–XI (Fig. 2), respectively. In the class VIII experiment we obtained cells with axial symmetry, an important feature since such cells will form tissues by the composition of layers of cells (i.e. the skin epithelium). About a fifty percent of the cells have a thin membrane, exhibiting the remainder fifty percent a large membrane. Most of the evolved cells have a cell nucleus. The class IX is composed by evolved cells resembling neurons (Fig. 4b). One of the main features of neuronal cells is they may have ramified structures or dendrites in order to enhance their connection via synapse with other neurons, receiving, processing and transmitting signals. The evolved neuronal biomorphs exhibited a thin and ramified membrane, bearing some of them a strong resemblance with the soma or perikaryon (the cell body) of real neurons. Similarly, some interesting results were obtained in class X. In the class X experiment, biomorphs resembling intestinal cells were evolved, but only during seven generations since the generation time was too long and the obtained forms were stables. Such cells form the intestinal tissue with folded structures in their membrane promoting the nutrients absorption. Most of the organelles are located at the top section of the cells. In the experiment the evolved cells were biomorphs with a thin membrane presenting folded structures and most organelles located near to such structures and not in the center of cells. Finally, erythrocytes or red blood cells were obtained under class XI experiment. Erythrocytes are specialized in the transport of oxygen to the other tissues, evolving in the present experiment rounded biomorph cells without cell nucleus (red) like real erythrocytes.

5. Discussion In this paper we have examined how unicellular organisms, as well as some specialized eukaryotic cells, evolve through a master plan of transformations under Darwinian natural selection from an initial population of cells simulated as Pickover biomorphs (Figs. 1 and 2). An interesting finding of our study is that the Pickover criterion promotes the evolutionary appearance of cell populations with a higher diversity of size and form than those populations of cells evolved under the Euclidean norm, the usual criterion applied to Julia sets. One of the basic principles of this paper is that we are mapping from a complex variable function to a set of morphologies. Taking into consideration our simulation results and those obtained by (Levin, 1994b), the biomorph’s form depends on the initial point Z0 as well as the complex function terms. Therefore, the complex function terms play different roles, particularly the complex variable and complex constant. For instance, Zn variable drives the population to rounded cell forms, whereas the trigonometric terms,

28

N.S. Mojica et al. / BioSystems 98 (2009) 19–30

i.e. sin(z), cos(z), and exponential term ez convey the population to horizontal cells or tubular forms. Further, we have found an important difference in the way that polynomials and trigonometric/exponential functions fill the cells with ‘organelles’, in particular through tessellation and parallel layers, respectively. In agreement with Levin (1994b), an important observation is the role of the complex constant c, showing a general tendency to disrupt the symmetry. Many biomorphs exhibit radial symmetry, others bilateral symmetry, etc. It is worth to note that in our model crossover or recombination is only carried out between those chromosome segments codifying for the value of the complex constant c. According to this tentative explanation the source of variability in the model should be the values of the real part (p) and imaginary part (q) of the complex constant c (or p + qi). Moreover, c could be modeling the expression or regulation of genes governing symmetry, e.g. Hox, Cnox, etc. (Ball et al., 2007; Schierwater et al., 2009). Nevertheless, the biomorph’s form is not only explained by the initial point Z0 and c, it is also (and mainly) explained by the convergence criterion (Jovanic et al., 2009). For instance, the cilia and appendages seen in biomorphs are absent when a relaxed convergence criterion |Re(z)| < ε or |Im(z)| < ε is substituted by a standard 

one, i.e. Re(z)2 + Im(z)2 < ε. In fact, the study of the relationship between the convergence criterion and the existing organismic morphologies was initiated by Pickover (1986, 1990, 1992) introducing the epsilon cross technique (Carlson, 1996). This method reveals the internal and external dendritic structures of a fractal as we conducted in the present paper. Is worthy of note, when convergence criterion is a root finding method (Levin, 1994b; Pickover, 1989), e.g. Newton and Halley methods, the obtained images do not bear a resemblance with organismic morphologies. In consequence, biomorph’s form is an ‘orbit trap’ or as we name in this paper an Evolutionary Stable Morphology (ESM). Thus, ‘organisms form’ is an attractor, but the ‘phenotype’ (organismic form and form features) that make up the attractor come about under a relaxed convergence criterion. This notion of biological form as attractor is similar to the “plant” attractor or “neural cell” attractor examples introduced by Kocic (2001) applying a tool named AIFS (Affine Invariant Iterated Function System). On the assumption that biomorphs’ form depends on the features addressed above, another important concern is about the robustness of the model. In agreement with Levin (1994b), when parameters are slowly changed in a polynomial function, e.g. from n = 2 to 10 in a Zn function, then the resulting morphology is always a circular form with spikes. In this study the author concludes that morphologies form “natural kinds” which could be equivalent to the cellular classes we have defined in this paper. Actually, what we have simulated with a genetic algorithm is the evolution of a morphogenetic field yielding from an undifferentiated unicellular form to diversified unicellular organisms. Consequently the morphogenetic field evolution is equivalent to a Darwinian ‘bauplan change’. One of the limitations of this study is the fact that the model and computer simulations cannot take into account the recombination in other genes than those four genes codifying for the value of the complex constant c. Another limitation is given by the way that a complex function P(z) is represented in the genetic algorithm. In this paper the complex function was codified as a chromosome instead of a Lisp tree, which is a more flexible representation as usually done in genetic programming (Koza, 1992). Also, the different values of the recombination and mutation probabilities were no tested, since in other related papers (Perales-Graván and Lahoz-Beltra, 2004, 2008) we had found that such values have no effect on the population evolution when genes have real values or represent a pattern (i.e. biomorphs, zebra stripes, etc.). It is interesting to observe a fast evolutive convergence to the optimum

cell in each one of the classes since, no matter the type of cell, the optimum is reached in a few generations. In our model, crossover or recombination is only carried out in the last four genes of the chromosomes. Since these genes codify for the value of the complex constant c, such constant would play an important biological role during morphogenesis (as we have already suggested, c could be modeling the expression of genes governing symmetry). This is congruent with Levin’s (1994a), where he states that c would represent some existing identical information that all of the cells share when they begin morphogenesis. This author also interprets that the obtained images, the biomorphs, are pre-pattern, thus bearing a resemblance with experimental micrographs depicting gene expression (obtained using a method that is known as in situ mRNA hybridization). Nevertheless, our opinion is that we obtain images of biomorphs that represent the final pattern of cells. In any case, many questions are still open. For instance, what would happen if |Re(z)| < ε AND |Im(z)| < ε would be the convergence criterion? Another interesting possibility would be to simulate a morphogenetic field A based on quaternion Julia fractals. In such a case we could obtain 3D biomorphs with the dimension reduction method described by Bourke (2001). As for the comparison with other theoretical models, what could serve as a common ground is the assumption of a 2D morphogenetic field. This field is defined by X-axis and Y-axis being the main difference among models the meaning of both axes. Thus, whereas Goodwin’s model (Goodwin, 1994; Goodwin and Trainor, 1985) is focused on the cortical cytoplasm, our model pays attention to the whole cytoplasm. In Goodwin’s model, the morphogenetic field is a mechano-chemical field, modeling the cytoskeleton–calcium coupling. In this approach, X-axis is the cytogel viscosity (including structural elements, e.g. microtubules) and Y-axis the kinetics of Ca2+ into the cytogel. Compared with our model we are focused on a mechano-fractal field, acting the Julia set and Pickover fractals as attractor of complex variable equations modeling ‘cellular bauplans’. Other models, i.e. Beloussov and Grabovsky (2007), assume that shape formation is conducted on a mechanical-stressed field formalized with a version of the classical Vander der Pol equations. That is, X-axis represents a slow mechanical stress whereas Y-axis a fast stress associated with Ca2+ dynamics and its effects on microtubules (Hotani et al., 1992). Curiously, Goodwin and Beloussov–Grabovsky proposals share their inspiration in marine life forms, the unicellular marine alga Acetabularia and the marine hydroid Dynamena, respectively. Even when both models have in common the assumption of a mechanical field and emphasize the role of calcium-microtubules interaction, the Beloussov–Grabovsky approach considers the notion of biomorphic forms (i.e. symmetric morphologies) as attractors. The idea of quasi-Julia set or Pickover attractors resembling biological morphologies or biomorphs is introduced by Pickover. However, is Levin who takes a step forwards connecting the Pickover approach with two relevant concepts in morphogenesis: the idea of Wolpert’s positional information and Turing’s diffusion–reaction model. Thus, in Levine model the morphogenetic field is a positional field modeled as a Pickover biomorph, with X and Y axes representing chemicals or morphogen concentrations. The complex variable equations are representing in Levin’s model the interaction between X and Y genes, expressing morphogenes. In contrast, eventhough our model shares some basic formal aspects with Levine model, it introduces “evolvability” in its formal structure. Of course, computer simulation experiments evolving Julia sets with genetic algorithms (Ashlock and Jamieson, 2007; Goertzel, 1995) are not a novelty, but in such experiments their authors don’t search for any biological issue relating fractals as Evolutionary Stable Morphologies (ESMs) with differentiation, growth and form of organisms. In summary, the current model describes how the cellular protoplasm and form of unicellular organisms and specialized

N.S. Mojica et al. / BioSystems 98 (2009) 19–30

eukaryotic cells can be evolved as Pickover biomorphs using a simple genetic algorithm. It also supports, under a theoretical realm, the existence of some attractors representing the functional and stable form of a cell—as stated, genuine cellular bauplans. In consequence, the morphogenetic field is attracted towards one or another attractor depending on the environment, modeling different environments with the fitness function. Additionally, with a refinement of the model we could attempt to simulate cancerous cells representing pathological states in the transformations of a cellular organism, as well as the normal development of animal and plant tissues; the plausibility of incorporating this distinction is an important point in favor of future developments of the model. Although we have addressed the problem of morphogenesis based on a positional information paradigm, in fact the simulation results also promote the views of Waddington, Goodwin and D’Arcy Thompson, which consider organisms as dynamical systems that evolve through a ‘master plan’ of transformations, and the classical notion of bauplan has been extended down to the cellular level. Theoretically, it is interesting that by means of standard mathematico-computational tools we have interconnected central ideas belonging to widely separated schools of biological thought. The proposed model elegantly connects with very recent ideas on mechanical morphogenetic field (Ainsworth, 2008). The classical view assumed a ‘chemical’ environment or field where signalling proteins and chemical gradients exclusively produce cell development and differentiation. However, the new empirical appreciations are throwing new light on the relevance of a cell ‘mechanical’ environment. In fact, experimental and theoretical papers had already contemplated the plausibility of a mechanical morphogenetic field. For instance, physical processes (cell wall expansion, diffusion) via cytoskeleton and with the intervention of other mechanisms (morphogens or electric fields) are though to be major aspect influencing the plant morphogenesis (Dupuy et al., 2006) and pollen aperture pattern (Ressayre et al., 1998). Studied gastropod mollusc eggs reveal a morphogenetic field related with the underlying cytoskeleton (Tyler et al., 1998). Lung morphogenesis has been modelled (Nelson and Manchester, 1988) using fractal geometries assuming that lung structure is largely determined by mechanical forces. It is also interesting to note how reaction–diffusion (morphogens) theory and its simulation have been coupled to growth explained by mechanical phenomena (Harrison et al., 2001). Further, the organization of the axial structure of amphibian embryo during morphogenesis is suggested (Beloussov, 1980) to be maintained by ‘tensile fields’ created by polarized cells. Also in vertebrates Gordon (2001, 2006) found physical ‘differentiation waves’ (contraction and expansion waves) traversing embryo; the most spectacular of these waves being the ectoderm contraction wave, leaving in its wake the neural plate. Notwithstanding that, the cytoskeleton is more than what the mathematical idea of the mechanical force field may imply. As Ainsworth has intuited (2008) it is also a locus where the informational transition between the local and the global may be achieved thanks to the “tensegrity” property. The way local chemical signals may be translated to mechanical asymmetries, and be broadcasted globally up to cell (and tissue) dimension, with a correlated quasiinstantaneous change of shape, represents a new way of thinking on the processes of development and physiology. It is here where the connection with new theoretical approaches to the “abstract” nature of cellular form is needed. We think that the cellular bauplan model, herein proposed, which is based on Julia sets and Pickover biomorphs, becomes the natural avenue to connect with the developmental field notion and the mechanobiological approach. At least, this is an attempt which will be continued by the authors in future works.

29

Acknowledgments ˜ We would like to thank I. Barradas Bribiesca and A. Cesena ˜ Quinones for useful comments and suggestions. The first author was supported by the University of Guanajuato (Mexico), the University of Guerrero (Mexico), CONACYT (Consejo Nacional de Ciencia y Tecnología, Mexico) as well as the Complutense University of Madrid (Spain). The second and third authors acknowledge support of the Instituto Aragonés de Ciencias de la Salud (Gobierno de Aragón, Spain); and the forth author was supported by Laboratorio de Bioinformatica, Complutense University of Madrid. Thanks to the anonymous reviewers for their suggestions and improvements.

References Ainsworth, C., 2008. Stretching the imagination. Nature 456, 696–699. Ashlock, D., Jamieson, B., 2007. Evolutionary exploration of generalized Julia sets. In: Proceedings of the 2007 IEEE Symposium on Computational Intelligence in Image and Signal Processing (CIISP 2007), pp. 163–170. Ball, E.E., M de Jong, D., Schierwater, B., Shinzato, C., Hayward, D.C., Miller, D.J., 2007. Implications of cnidarian gene expression patterns for the origins of bilaterality—is the glass half full or half empty? Integr. Comp. Biol. 47, 701–711. Basler, K., Struhl, G., 1994. Compartment boundaries and the control of Drosophila limb pattern by hedge hog protein. Nature 368, 208–214. Beloussov, L.V., 1980. The role of tensile fields and contact cell polarization in the morphogenesis of amphibian axial rudiments. Roux’s Arch. Dev. Biol. 188, 1–7. Beloussov, L.V., Grabovsky, V.I., 2007. Information about a form (on the dynamic laws of morphogenesis). BioSystems 87, 204–214. Blickle, T., Thiele, L., 1995. A mathematical analysis of tournament selection. In: Proceedings of the 6th International Conference on Genetic Algorithms, Morgan Kaufmann, Pittsburgh, PA, pp. 9–16. Bourke, P., 2001. Quaternion Julia fractals. http://local.wasp.uwa.edu.au/∼pbourke/ fractals/quatjulia/. Briere, C., Goodwin, B.C., 1990. Effects of calcium input/output on the stability of a system for calcium-regulated viscoelastic strain fields. J. Math. Biol. 28 (5), 585–593. Carlson, P.W., 1996. Pseudo-3-D rendering methods for fractals in the complex plane. Comput. Graph. 20 (5), 751–758. Chovnjuk, Y.V., Ovsyannikova, T.N., Ivanovskaya, A.V., 2001. Stochastic branching process model of phylogeny in mm wave range electromagnetic field. In: The Fourth International Kharkov Symposium on Physics and Engineering of Millimeter and Sub-Millimeter Waves, IEEE, vol. 2, Kharkov, Ukraine, pp. 949–951. Davis, L. (Ed.), 1991. Handbook of Genetic Algorithms. Van Nostrand Reinhold, New York. Dawes-Hoang, R.E., Parmar, K.M., Christiansen, A.E., Phelps, C.B., Brand, A.H., Wieschaus, E.F., 2005. Folded gastrulation, cell shape change and the control of myosin localization. Development 132 (18), 4165–4178. De Pedro, M.A., 2004. Topological domains in the cell wall of Escherichia coli. In: Vicente, M., Tamames, J., Valencia, A. (Eds.), Molecules in Time and Space: Bacterial Shape, Division and Phylogeny. Springer, New York, pp. 27–58. Devaney, R.L., 1994. The complex dynamics of quadratic polynomials. In: Complex Dynamical Systems: The Mathematics Behind the Mandelbrot and Julia Sets. Oxford University Press, Oxford, pp. 1–27. Dupuy, L., Mackenzie, J.P., Haseloff, J.P., 2006. A biomechanical model for the study of plant morphogenesis: Coleocheate orbicularis, a 2D study species. In: 5th Plant Biomechanics Conference, Stockholm. Frankel, J., 1975. Pattern formation in ciliary organelle systems of ciliated protozoa. In: Cell Patterning. Ciba Fdn Symp., Elsevier, London, pp. 25–49. Gierer, A., 1981. Generation of biological patterns and form: some physical, mathematical, and logical aspects. Progr. Biophys. Mol. Biol. 37, 1–47. Gierer, A., Meinhardt, H., 1972. A theory of biological pattern formation. Kybernetik 12, 30–39. Goertzel, B., 1995. Rapid generation of strange attractors with the eugenic genetic algorithm. Comput. Graph. 19 (1), 151–156. Goldberg, D.E., 1989. Genetic Algorithms in Search, Optimization and Machine Learning. Addison-Wesley, Reading, MA. Gordon, R., 2001. Making waves: the paradigms of developmental biology and their impact on artificial life and embryonics. Cybern. Syst.: Int. J. 32, 443–458. Gordon, R., 2006. Mechanics in embryogenesis and embryonics: prime mover or epiphenomenon? Int. J. Dev. Biol. 50, 245–253. Goodwin, B.C., Trainor, L.E.H., 1985. Tip and whorl morphogenesis in Acetabularia by calcium-regulated strain fields. J. Theor. Biol. 117, 79–106. Goodwin, B.C., 1994. How the Leopard Changed its Spots: The Evolution of Complexity. Scribner, New York. Gurwitsch, A., 1922. Über den Begriff des embryonalen Feldes. Roux Arch. Ent. Org. 51, 383–415. Harold, F.M., 1990. To shape a cell: an inquiry into the causes of morphogenesis of microorganisms. Microbiol. Rev. 54, 381–431. Harrison, L.G., Wehner, S., Holloway, D.M., 2001. Complex morphogenesis of surfaces: theory and experiment on coupling of reaction–diffusion patterning to growth. Faraday Discuss 120, 277–294.

30

N.S. Mojica et al. / BioSystems 98 (2009) 19–30

Hotani, H., Lahoz-Beltra, R., Combs, B., Hameroff, S., Rasmussen, S., 1992. Microtubule dynamics, liposomes and artificial cells: in vitro observation and cellular automata simulation of microtubule assembly/disassembly and membrane morphogenesis. Nanobiology 1 (1), 61–74. Jovanic, R., Tuba, M., Simian, D., 2009. Recent Advances in Mathematics and Computers in Biology and Chemistry. Proceedings of the 10th WSEAS International Conference on Mathematics and Computers in Biology and Chemistry (MCBC’09), Prague, Czech Republic, pp. 162–167. Kavlock, R.J., Setzer, R.W., 1996. The road to embryologically based dose–response models. Environ. Health Perspect. 104 (1), 107–121. Kocic, L.M., 2001. AIFS—a tool for biomorphic fractal modeling. Nonlinear Dynam. Psychol. Life Sci. 5 (1), 45–63. Koza, J.R., 1992. Genetic Programming: On the Programming of Computers by Means of Natural Selection. The MIT Press, Cambridge, MA. Lahoz-Beltra, R., 2001. Evolving hardware as model of enzyme evolution. BioSystems 61, 15–25. Lahoz-Beltra, R., 2004. Bioinformatica: Simulación, Vida Artificial e Inteligencia Artificial. Ediciones Díaz de Santos, Madrid (Transl.: Spanish). Lahoz-Beltra, R., 2008. ¿Juega Darwin a los Dados? Editorial NIVOLA, Madrid (Transl.: Spanish). Lahoz-Beltra, R., Selem Mojica, N., Perales-Gravan, C., Navarro, J., Marijuán, P.C., 2008. Towards a Morphogenetic Field Theory. http://bioinformatica.net/mfproject/index.html. Levin, M., 1994a. A Julia set model of field-directed morphogenesis: developmental biology and artificial life. Comput. Appl. Biosci. 10 (2), 85–103. Levin, M., 1994b. Discontinuous and alternate Q-system fractals. Comput. Graph. 18 (6), 873–884. Limozin, L., Denet, B., Pelce, P., 1997. Ion currents generated by tip growing cells. Phys. Rev. Lett. 78 (25), 4881–4884. Lodish, H., Berk, A., Zipursky, S.L., Matsudaira, P., Baltimore, D., Darnell, J., 2000. Molecular Cell Biology. W.H. Freeman and Company, New York. Lutton, E., 1999. Genetic Algorithms and Fractals—Algorithmes génétiques et fractales, Spécialité Informatique, Habilitation à diriger des recherches. Université Paris XI Orsay, 11 Février. Lynn, D.H., Tucker, J.B., 1976. Cell size and proportional distance assessment during determination of organelle position in the cortex of the ciliate Tetrahymena. J. Cell Sci. 21 (1), 35–46. Marijuán, P.C., 1996. Information and symmetry in the biological and social realm: new avenues of inquiry. Symmetry: Cult. Sci. 7 (3), 281–294. Meinhardt, H., 2003. The Algorithmic Beauty of Sea Shells, 3rd enlarged edition. Springer (with PC—software), Heidelberg, New York. Nelson, T.T., Manchester, D.K., 1988. Modeling of lung morphogenesis using fractal geometries. IEEE Trans. Med. Imag. 7 (4), 321–327. Nuccitelli, R., 1984. The involvement of transcellular ion currents and electric fields in pattern formation. In: Pattern Formation—A Primer in Developmental Biology. Macmillan Publishing Company, London, pp. 23–46. Nuccitelli, R., 1988. Ionic currents in morphogenesis. Experientia 44, 657– 666. O’Shea, P.S., 1988. Physical fields and cellular organization: field-dependent mechanisms of morphogenesis. Experientia 44, 684–694.

Perales-Graván, C.P., Lahoz-Beltra, R., 2004. Evolving morphogenetic fields in the zebra skin pattern based on Turing’s morphogen hypothesis. Int. J. Appl. Math. Comput. Sci. 14 (3), 351–361. Perales-Graván, C.P., Lahoz-Beltra, R., 2008. An AM radio receiver designed with a genetic algorithm based on a bacterial conjugation genetic operator. IEEE Trans. Evol. Comput. 12 (2), 129–142. Perry, J.J., Staley, J.T., Lory, S., 2002. Microbial Life. Sinauer Associates, Publishers, Sunderland, MA. Pickover, C.A., 1986. Biomorphs: computer displays of biological forms generated from mathematical feedback loops. Comput. Graph. Forum 5, 313–316. Pickover, C.A., 1989. Chaotic fragmentation in Halley’s paradise. Physica Scripta 39, 193–195. Pickover, C.A., 1990. Computers, Pattern, Chaos and Beauty. St. Martin’s Press, New York. Pickover, C.A., 1992. Computers and Imagination. St. Martin’s Press, New York. Raff, R.A., 1996. The Shape of Life: Genes, Development, and the Evolution of Animal Form. The University of Chicago Press, Chicago. Ressayre, A., Godelle, B., Mignot, A., Gouyon, P.H., 1998. A morphogenetic model accounting for pollen aperture pattern in flowering plants. J. Theor. Biol. 193, 321–334. Schierwater, B., Eitel, M., Jakob, W., Osigus, H.-J., Hadrys, H., Dellaporta, S.L., Kolokotronis, S.-O., DeSalle, R., 2009. Concatenated analysis sheds light on early metazoan evolution and fuels a modern “Urmetazoon” hypothesis. PLoS Biol. 7 (1), 36–43. Spemann, H., 1921. Die Erzeugung tierischer Chimaeren durch heteroplastiche embryonale Transplantation zwischen Triton cristatus u. taeniatus. Wilhelm Roux Arch. Entwicklungsmech. Org. 48, 533–570. Spemann, H., Mangold, H., 1924. Über Induktion von Embryonanlagen durch Implantation artfremder Organisatoren. Roux’ Arch. f. Entw. mech. 100, 599–638. Spemann, H., 1938. Embryonic Development and Induction. Yale University Press, New Haven. Thompson, D’Arcy W., Bonner, J.T., 1992. On Growth and Form: The Complete Revised Edition. Cambridge University Press, Cambridge. Turing, A.M., 1952. The chemical basis of morphogenesis. Philos. Trans. R. Soc. Lond. Ser. B, Biol. Sci. 237, 37–72. Tyler, S.E.B., Butler, R.D., Kimber, S.J., 1998. Morphological evidence for a morphogenetic field in gastropod mollusc eggs. Int. J. Dev. Biol. 42, 79–85. Vences, L., Rudomin, I., 1997. Genetic algorithms for fractal image and image sequence compression. In: Proceedings Computación Visual, Universidad Nacional Autónoma de Mexico, pp. 35–44. Waddington, C.H., 1940. Organisers and Genes. Cambridge University Press, Cambridge. Waddington, C.H., 1942. The epigenotype. Endeavour 1, 18–20. Waddington, C.H., 1957. The Strategy of Genes. Allen & Unwin, London. Weiss, P., 1923. Naturwissen 11, 669. Weiss, P., 1939. Principles of Development. Holt and Co., New York. Wolpert, L., 1969. Positional information and the spatial pattern of cellular differentiation. J. Theor. Biol. 25, 1–47. Young, D.A., 1984. A local activator-inhibitor model of vertebrate skin patterns. Math. Biosci. 72, 51–58.