extfieldsofstudy:
Computer Science
Medicine
Biology

text:
Origins of the Embryo: Self-organization through cybernetic regulation.

The construction of an embryo from a single cell precursor is a highly complex process. Evolutionary emergence of the first embryos is even more complex, and involves both a transition to multicellularity along with the establishment of developmental mechanisms. We propose that embryogenesis relies on a community of cells conforming to a regulatory model of emergent multicellularity. This model draws together multiple threads in the scientific literature, from complexity theory to cybernetics, and from thermodynamic entropy to artificial life. All of these strands come together to inform a model of goal-oriented regulation for emergent structures in early life. This is an important step in the evolution of early life, as well as the emergence of complex life in the earliest habitats. Our model, called the cybernetic embryo, allows for a systems-level view of the embryogenetic process.

The Origin of a Cybernetic Embryo
"Deep understanding of causality sometimes requires the understanding of very large patterns and their abstract relationships and interactions, not just the understanding of microscopic objects interacting in microscopic time intervals." Douglas Hofstadter, I am a Strange Loop (Hofstadter, 2007).
The idea of applying cybernetics to understand embryogenesis has been proposed as an alternative means of understanding the processes of self-reference, collective organization, and differentiation during tissue formation and organogenesis (Gordon and Stone, 2017). According to this view, a series of nested cybernetic (self-regulating) systems are responsible for generating differentiated tissues and multicellular structures. This model exists within a larger, more complex picture of differentiation, which includes experimental work in frogs and sea urchins (Kearl, 2012), theoretical work by Gordon (1999), and even more recent work on electro-chemical signaling networks (Levin and Pezzulo, 2016). The nested structures that result from a cybernetic developmental process includes both hierarchical and heterarchical configurations (Bruni and Giorgi, 2015;von Goldammer et al., 2003). While many biological systems are hierarchical, heterarchies also occur in contexts where the constituent elements of a system can be configured in a number of equivalent ways (McCullough, 1945). These complex systems also utilize signals from outside the embryo, which results in an open thermodynamic system. How does cybernetics describe highly complex biological systems? Cybernetics is the study of the predictive behavior of systems as they converge upon goal states. As a general theory of systems, cybernetics provides an abstract and holistic view of a biological system's universal characteristics. It also provides a simplified dynamical description of these biological systems as 1 they transition between equilibrium and non-equilibrium states (Weiner, 1961;Ashby, 1963;Pask, 1961). A cybernetic view of embryogenesis suggests an inviolate level at both the top and the bottom of the analysis: at the bottom the origin of organismal complexity, at the top the development of the embryo in successive stages in evolution. Ultimately, the central theme of the cybernetic embryo model is in providing a formal description of how regulatory mechanisms might behave in the regime between order and chaos (Gelens et al., 2015;McIsaac et al., 2011).
In this chapter, we will look beyond the specificity of the biological processes related to processes such as differentiation or induction. Rather, we will ask why groups of totipotent cells that later differentiate in a coordinated fashion are able to produce high degrees of complexity. We will approach this by treating the embryo as an abstract set of processes that represent the deep structure of developmental processes. We will establish a development model based on insights from cybernetics, game theory, computational linguistics, biological information, and physics that can explain the process of developmental history at both the cellular and organismal levels. A developmental model based on the principles of cybernetic regulation could become quite unwieldy. Therefore, we must state the requirements for origination of the cybernetic embryo and its processes. This leads us to two interrelated questions: 1) how do multicellular systems emerge, particularly in a developmental/life-history context, and 2) what is the fundamental unit of regulation in the embryo?
One primary concern of cybernetic embryo origins is the origin of regulative processes themselves, particularly as they relate to the origin of multicellular systems. Broadly speaking, the embryo satisfies two criteria in defining living systems: embryos are self-sustaining (Luisi, 2003) and self-replicating (Schuster and Sigmund, 1983). Whether found within eggs or within seeds, embryos require a supportive environment. In this sense, embryos are not completely selfsustaining. However, Turchin's (1995) definition of metasystems suggests that an embryogenesis embedded in the proper physiochemical environment can create and support the embryo as a selfsustaining entity. The criterion of self-sustenance is also related to the process of autopoiesis (Varela et al., 1974). In Abel (2009), autopoiesis as defined in the embryo consists of three parts: a network of reactions (or interactions) capable of regenerating all of a system's components existing within a self-defined boundary. In the embryo, we have both totipotent and pluripotent cells as well as a directed form of regeneration in which the basic arrangement (e.g. organismal structure) of the parent is recapitulated in the offspring. This recapitulation occurs through a set of feedback mechanisms (Varela et al., 1974) between genome, organism, and environment.

Mechanisms of the Cybernetic Embryo
Two ideas that will run through this chapter are the Law of Requisite Variety and the Good Regulator Theorem. Taken collectively, these ideas will provide a link between the cybernetic embryo as a conceptual model and the dynamic aspects of evolution and development. The Law of Requisite Variety, or the first law of cybernetics, defines the role of variation in a regulated system (Ashby, 1963). In an equilibrium system, variation is advantageous in terms of both competition and controllability. Particularly in terms of controllability, variation of the controller must be greater than the number of states controlled (Ashby, 1963). In this way, cybernetic systems can exhibit great robustness without having to predict every possible state in the system. While this is similar to the anticipatory systems of Rosen (1985), we have integrated such an argument with the notion of epigenetic landscapes (Waddington, 1935;Waddington, 1962) to yield a regulatory model that is both dynamic and generative.
In the case of our cybernetic embryo, we propose that the determinant of all developmental trajectories (evolution) produces an infinitely greater variety of possible states than the developmental process can encompass in any one species. In some cases, such as in the nematode Caenorhabditis elegans, development is mostly deterministic in terms of pathways taken during the developmental process (Penigault and Felix, 2011). In this case, evolution leaves no room for developmental contingencies, but produces a highly controlled trajectory. In other cases, such as the phenomenon of phenocopying (phenotypic diversity from the same genotype) in mammalian development (Sharma et al., 2016), evolution provides an avenue for robustness by leaving a number of possible contingencies for development that result in the same phenotype. From a cybernetics perspective, evolution has now resulted in a finite set of strategies (Koestler, 1967) expressed in the form of unique developmental trajectories. However, the amount of variation available to the evolutionary process influences what control strategies are available with relation to the developmental system. If the genetic pathway controlling development is highly complex, then we should expect more diversity in the developmental process. Overall, controllability (for definition in complex networks, see Liu et al., 2011) will vary based on the amount of organismal diversity provided by evolution and the number of viable developmental pathways.
Understanding this relationship in terms of the biology of an embryo also requires an appreciation for the Every Good Regulator Theorem (EGRT; see Conant and Ashby, 1970), which states that well-regulated systems must be a model of the system. In other words, regulatory mechanisms are meta-structural, not necessarily linked to specific biological components. We propose that in the embryo, the model of regulation is a discrete combinatorial set of pathways that connect the products of evolution and the instructions of development.

Candidate Processes
In a wide range of biological cell types, and even in protocells (Chen and Walde, 2010), we can see that single cells possess many regulatory mechanisms. In fact, the term "regulator" might seem too general given the complexity of cell function. Yet in this case, the term "regulator" is meant to be more general, as in a mechanism that allows for cell lines to replicate and remain self-contained as a living system. While our concern is largely theoretical, it is also related to the transition from a single cell to multiple cells in the same organismal system. In this sense then, we are also interested in how clusters of cells are able to regulate themselves as a single self-contained entity. When control occurs at the global level, it influences all other function in the cellular unit. Depending on the context, the regulatory mechanism can change scope during the process, with the same mode of control occurring at different scales (Lodhi and Muggleton, 2010). We will discuss the possible mechanisms for control later in the chapter. For our present concern, it is important to keep the black box regulator distinct from the theoretical construct.
Related to this is the need to understand what makes for an efficient regulator. Many processes in the cell are efficient, yet in terms of the theoretical regulator, the measure for efficiency is quite different. This is associated with the scale of regulation. Let us use a specific metabolic pathway as an example of a well-defined (non-theoretical) regulator. While this pathway can be considered an efficient regulator, its scale is not relevant to the entirety of the organismal (single or multicell) system. It is the singular course-grained regulator at the macroscale, which must be efficient for our cybernetic embryo to be persistent and self-sustaining.
Upon presenting the case for effective biological regulation, we can now consider the spontaneous order of embryos. More specifically, we can consider emergent order as it relates to the first embryos, which themselves are examples of a single cell to multicell transition in lifehistory (Kalinka and Tomancak, 2012). In a sense, embryogenesis encodes the evolutionary transition of multicellular complexity into a formal biological process. How exactly does the unfolding of this process become a quasi-repeatable phenomenon? In short, embryogenesis proceeds through a process of instruction resulting from genomic content and the spatial positioning of individual cells. The process is quasi-repeatable because while embryogenesis often yields similar results in embryos of the same species, depending on the species and mode of development there is some room for error and variation. Nevertheless, the answer we are looking for already exists in the epigenetic landscapes of Conrad Waddington (Waddington, 1957;Morris et al., 2016;Goldberg et al., 2007).
The epigenetic landscape is a tree-like structure that favors certain developmental pathways over others. Waddington referred to these favored pathways as canals (Siegal and Bergman, 2002), which when reinforced across individuals and over generations leads to the canalization of the developmental process (Waddington, 1957). The landscape is a discrete possibility space (based on the geometry of a Hamiltonian phase space) representing an exhaustive number of possible developmental outcomes (Goldberg et al., 2007). While Waddington originally meant for the epigenetic landscape to be metaphoric and descriptive in nature, more recent models envision development as a process of self-assembly, driven by free energy minimization and the ability of cells to reproduce a target shape in a collective manner (Friston et al., 2015).
It is our contention that both historical contingency and energetic constraints are responsible for determining the shape of canalization through this landscape. Specific developmental events necessarily lead to a finite number of subsequent events. This presents itself as a form of epigenetic memory, which when combined with specific instructions, can produce a tightly constrained tree topology along this landscape. We can observe this process in the developmental founder cells of Caenorhabditis elegans (Rose and Kemphues, 1998) or the establishment of developmental modules in vertebrates (Kuratani, 1998). The epigenetic landscape construct is also useful for understanding contexts in which these developmental constraints no longer act to focus the multicellular system. For example, we can examine the difference between increasing degrees of freedom and reduction of degrees of freedom due to a one-way embryogenetic development (for general concept, see Wissner-Gross and Freer, 2013). In a strictly developmental context, the available degrees of freedom are reduced as the embryo becomes increasingly specified. However, for biological processes such as cancer or aging, the degrees of freedom are increased.
Intuitively, the addition of cells to an organism might suggest a sort of combinatorial explosion of complexity that is biologically intractable. Epigenetic constraints resulting from developmental contingency (Oyama et al., 2003) suggest that the opposite holds true. In the embryo, the establishment of contingency results in part from the hierarchical organization of multicellular complexity and thermodynamic constraints. As a quasi-historical process, embryogenesis is shaped by developmental contingencies, which limit the role of environmental variation and are enforced by the twin phenomena of thermodynamic entropy and genetic information. As part of a more abstract cybernetic process, contingencies serve to shape the unfolding of information itself over time. As closed-ended dynamic systems (such as embryogenesis) progress, messages are decoded and assembled within an increasingly specific context (Weiner, 1954).

Thermodynamics and Complexity
In the middle of 19 th century, the irreversible process was first identified as a subject of theoretical interest, and defined a new subdiscipline of physics called thermodynamics. Later, Rudolf Clausius characterized the Second Law of Thermodynamics as the inevitable growth of entropy in any closed system. Ludwig Boltzmann developed a mathematical basis for the Second Law several decades later (Sharp and Matschinsky, 2015). In particular, Boltzmann's H-theorem that allows us to characterize thermodynamic irreversibility using a molecular-kinetic interpretation of heat processes and a comparison of thermodynamic state probabilities. A halfcentury later, a means to characterize the role of entropy in information was developed by Shannon (1948) and gave birth to a new scientific discipline revolving around the theory of information. Ben-Naim (2015) has remarked thusly on the relationship between the two: "Entropy, as it is now recognized, does not mean 'transformation', or 'change', or 'turn'. It does mean information". While many people accept the congruities between thermodynamic entropy and Shannon information, treating them as equivalent is not appropriate for all types of systems. In the concept of negentropy (Brillouin, 1953), entropy and information are treated as separate entities. According to this view, there is a subtle but important distinction made between entropy and information in cases where entropy works against the growth of complexity. In living systems, negentropy might act to limit the role of entropy but not information during the accumulation of complexity.
One way in which biological complexity might interact with information is by characterizing the information use in the regulation of processes as a form of algorithmic (or descriptive) complexity (Solomonoff, 1964;Kolmogorov;1965;and Chaitin, 1969). Algorithmic complexity is relatable mathematically to Shannon Information, and reveals direct dependency between a given system's complexity and its total degrees of freedom. Although not equivalent to Shannon information, algorithmic complexity and entropy are mathematically similar. These mathematical interrelationships allow us to generalize the Second Law of Thermodynamic as a general law of complification (Mikhailovsky and Levich, 2015). Taking into account conditional complexity derived through degrees of freedom as indicated above, this law might also be the way in which living systems maximize their freedom.
We see countless examples of complification in the world around us. Complification is seen in a universe developing from simplicity of Big Bang and baryon soup to its current complexity of metagalaxies, galaxies, stars, and planets. It is matter born of quarks, electrons, neutrinos, gluons, etc., and has come up to macromolecules, cells, and living organisms. Perhaps the best example of complification is the evolution of complexity encompassing prokaryotes, eukaryotes, and multicellular organisms. It is characteristic of each organism that develops from fertilized zygote to mature organism. It defines the evolution of humanity, from the simplicity of scattered tribes to complexity of the contemporary world with nations, corporations, and stock markets. These examples are far from being exhaustive. However, these examples clearly contradict the interpretation of complexity as defined by maximum entropic states and randomization (chaos) as simple processes. In fact, these are all complex processes for which the maximum degree of complexity lies on the edge of order and chaos (Lineweaver, 2013). We also see this with the growth of real-world algorithmic complexity in structures that age and/or evolve over time (Papadimitriou, 2014). Yet increases in algorithmic complexity usually do not reach their maximum values regardless of system.
The reason for this apparent suboptimality is simple: growth in complexity is limited by natural law. Laws of nature lead to a notion of permanent counteraction between the general law of complification and all other laws of nature. These laws are the only entities that do not obey the general law of complification and have not changed, as far as we know, during all the existence of our universe. This counteraction reveals itself at all the levels of hierarchy. Even the maximum entropy of ideal gas considered by Boltzmann in his H-theorem is not a real chaos. Its material points move in straight lines (Newton's First Law), collide elastically (Newton's Third Law), and their average energy, i.e., temperature, is proportional to the square of their velocity (Newton's Second Law). With increasing complexity and diversity, limitations by various laws of nature increase as well (Fuentes, 2014).
However, there is at least one way out that allows complexity to break through local barriers set by the laws of nature. This is called hierarchogenesis, and specific examples of hierarchogenetic events are described in Mikhailovsky and Gordon (2017). In the context of a cybernetic embryo, hierarchogenesis proceeds by systems of a given hierarchical level uniting from individual components. The total information of such a new super-system at a new level of hierarchy will be equal to the sum of information values for each of united systems (due to the additivity of entropy and information) plus information produced by interactions among these systems. This allows systems and their associated complexity to increase to new plateaus over time. While hierarchogenesis provides an emergent basis for the origin of such structures, increasing complexity may also be due to the law of requisite variety and emergent cybernetic control. In the next section, we will consider how this might be so.

Information and life
To understand how hierarchogenesis leads to actual biological complexity, we can turn to a process central to embryogenesis called morphogenesis. Morphogenesis is a process of diversifying organs and tissues within an organism to complete its life cycle. During the course of morphogenesis, hierarchical structures are created through processes of regulation. Yet the biology of phenotype and molecules demonstrate these regulatory relationships in an indirect manner. In fact, nature features many mechanisms, chemical, mechanical, electrical, etc., to accomplish this. At their core is the processes of information signaling and utilization.
To abstract this within a cellular automaton, the Morphozoic model (Portegys et al., 2017) has shown that tractable rule sets composed of hierarchically nested cellular neighborhoods are capable of producing both local and global information processing effects that can simulate reaction-diffusion morphogenesis (Turing, 1952), gastrulation, axon pathfinding, and other phenomena. Stephen Wolfram (2002) has demonstrated that simple rule sets can produce dynamic patterns within them, and that this complexity can be abstracted into hierarchical graphs (Wolfram, 2015).
We know that complex living systems lie in a middle ground between high and low entropy (Adami and Brown, 1994), and that perhaps regulatory mechanisms are responsible for keeping such systems within a critical set of values. As with the idea of negentropy, there is evidence that natural complex systems become more robust as they grow in complexity, and this robustness enables mechanisms such as self-sustenance and mutability. Even artificial systems seem to show this property under certain circumstances. For example, variations of circuits that can implement a logical function become more numerous as the number of components grows (Raman and Wagner, 2011). In fact, it may be that the evolvability of a complex biological system and its tendency for innovation are features of the information content itself (Ball, 2015;Eigen and Schuster, 1979).
Continuing along these lines, studies of genetic variations indicates that "fit" configurations tend to be located on particular "paths" in the space of all possible configurations (Kirshner and Gerhart, 2006). This allows one successful genotype to migrate through mutation into other successful ones that could be more complex, and can be taken as a refutation of the oftcited implausibility of life arising due to its sheer combinatorial complexity (Parter et al., 2008). Ball (2015) demonstrates this in discussing the interconnected nature of organismal complexity at multiple scales in determining fitness. For example, there are pathways in highly dimensional, complex systems that lead to "fit" phenotypes in only a few mutational steps. This is in part a property of combinatorics, where neural networks capable of sorting a very large number of inputs are constrained by a relatively small number of plausible pathways (Lin and Tegmark, 2016). In the biological case, however, greater complexity leads to more paths to fitness instead of making it increasingly less probable. In the context of the evolution of development, this provides a means to find stable paths to embryogenesis.

Emergence of Cybernetic Control via Evolutionary Transition and Constraint
In this section of our chapter, we will discuss how cybernetic elements can be put together to demonstrate the emergence of phenotypic complexity. We will mathematically and graphically define cybernetic elements, and then place such systems of elements into both a toy model and evolutionary context. Our toy model of choice for this demonstration is the Vehicle model of nervous system complexity, an approach introduced by Valentino Braitenberg (1984) to synthesize simple embodied nervous systems of varying complexities. To do this, we must represent the components of the vehicles as simple cybernetic elements. We must also consider what these cybernetic elements represent when arrayed in a set of interrelationships. By the end of this section, this will be explained and placed in the context of the cybernetic embryo.
While even the simplest cybernetic elements represent dynamic processes, they can also self-assemble into complex systems. Cybernetic systems can therefore represent changes in complexity of complex systems over time. This process of integration represents interactive processes ranging from symbiotic (Douglas, 2014) to cooperative (Pfeiffer and Bonhoeffer, 2003), and is representative of major evolutionary transitions (Maynard-Smith, 1983;Szathmary and Smith, 1995). In a less-abstract manner, this reintegration can represent the embryogenesis of a Vehicle, or transformation from a set of regulators functioning independently, to a set of autoregulating cells that work together to reproduce and function as higher-level units.
We can begin our discussion of cybernetic element complexity with simple cybernetic elements that can be represented both mathematically and graphically. We can characterize a set of cybernetic elements as an inferential structure, or a graphical model which models a specific hypothesis about regulation at a single biological scale. The simplest cybernetic structure is a firstorder feedback loop (Figure 1) over a unit of time, which can be defined with a set of partial differential equations where dt is the cybernetic structure over time, and ds is the cybernetic unit across space. The cybernetic structure consists of a unit x that receives a feedback loop (FB), a unit y that receives a feedforward loop (FF). Cybernetic units are used to characterize regulation independently over time (as is the case with serial cell divisions) or across space (as is the case in intercellular signaling).
The cybernetic models introduced here are adapted from models of molecular production and turnover (Alicea, 2014), and has been modified to approximate the regulation of discrete space-time processes. This differs from approaches that treat the embryo purely as a bag of chemical inducers and morphological gradients (Turing, 1952). In application to the cybernetic embryo concept, we combine simple cybernetic elements with pathways to serve as regulators of cellular state. Each set of cybernetic element and pathway combination (examples shown in Figures 1-3) represent idealized motifs, which simplifies the basic regulatory configuration independent of tissue type or even scale. With a focus on process rather than phenomenon, oversimplistic representations provide us with an idealized way to model data related to complex processes. Cybernetic elements as introduced here may or may not be implemented in a material fashion. As a purely theoretical tool, networks of cybernetic elements can provide structure to messy systems ranging from the dynamics of transcription in cells to the interactions between cells themselves. The use of cybernetic element as a fundamental unit is broadly applicable, and serve as a problem structuring method (Mingers and Rosenhead, 2004) to make complex biological systems tractable. However, black box cybernetic elements can also be represented as formal replicator systems such as von Neumann self-replicating automata (Sipper, 1998;von Neumann and Burks, 1966).
Therefore, we also propose that cybernetic elements can exist as ensembles, which in this context represent a collection of cells with transient or nonexistent functional interrelationships. Let us now assume that an ensemble of simple cybernetic elements can self-replicate at rate τ, but with a simple structure that does not interact with any of the other replicators. This results in a hybrid model of first-order regulation (serial) of a set of cybernetic elements with population dynamics. This can be defined using a form of the Quasispecies equation (Eigen and Schuster, 1979), and can be defined in the following manner where ̇ is a differential term, is the contribution of replicator to the production of replicator , is a specific replicant at a particular timepoint, is the replicator, and is a coefficient that represents the rate of degradation and/or mortality.
As both feedback in individual units and the size of the population fluctuates, it becomes obvious that systems of increased complexity will only emerge in a highly transient manner. Therefore, a higher level of regulatory control is required. This higher-order control need not be hierarchical, but does need to take advantage of the products resulting from spontaneous organization (Eigen and Schuster, 1979). This results from cooperative effects between cybernetic elements resulting from functional needs. One source of emergent interdependency involves the long-term consequences of regulation for each individual cell. Consider the time-series that result from Equation 1: in cases of either too much (runaway) feedback or too little (dampening) input, a single regulator can become dysregulated quite easily. In these cases, a more elaborate cybernetic element can form, and may result in a larger regulatory system (see Figure 2).
Larger regulatory systems generally result from the inclusion of unit along with corresponding feedforward and feedback connections. An example of this is analogous to the static nervous system design for Vehicle #2 in Braitenberg (1984). In this case (Figure 2), each feedback unit is cross inhibited, while each feedforward unit receives a first-order feedback from the crossinhibited feedback unit. This is mathematically defined as = ′⁄ , where dt is the cybernetic structure over time, ds is the cybernetic unit across space, x and x' are parallel feedforward units, y and y' are parallel feedback units, and E is the global systemic (e.g. behavioral) output.
Another form of interdependency resulting from regulatory needs results in the formation of futile (or metabolic) cycles (Qian and Beard, 2006). Futile cycles (see Figure 3) are indefinite first-order feedback loops that utilize two multiplicative noise channels: one in the feedforward loop, and one in the feedback loop (Samoilov et al., 2005). The feedforward noise channel is positive, while the feedback noise channel is negative. Such channels produce no tangible output, but are both energetically active and produce an entropic output. It is not clear whether futile cycles are adaptive, as they consume energy but produce no tangible output. As a part of larger systems, they serve a purely regulatory function (Tu and McKnight, 2006). This can be described mathematically as = x (σ) + y (−σ) = x (σ) + y (−σ)

[4]
where dt is the cybernetic structure over time, ds is the cybernetic unit across space, x(σ) is the noise-driven feedforward (FF) component, and y(-σ) is the noise-driven feedback (FB) component (which result in negative feedback from element y to element x). Futile cycles can evolve from first-order feedback loops as either a support structure for other, more complex regulatory elements, or as a transitive structure from first-order feedback loops to more complex elements. Thus, transitional elements such as the futile cycle may allow for the buildup of complexity amongst a set of cybernetic elements.
Biologically speaking, moving from a collection of simple cybernetic elements to a higherlevel system of interacting cybernetic elements is an evolutionary transition (Szathmary and Smith, 1995). While this could be concurrent with the transition from single-cell to multi-cell organisms, it could also characterize the transition between Prokaryotes and Eukaryotes or the transition from RNA world to genomes containing DNA. We can treat this variety of evolutionary transitions in such a generic fashion because the biological process can be demonstrated both graphically and quantitatively as the expansion of a generic regulatory system.
As cybernetic transitions allow for computational flexibility that scales to multiple problems and hypotheses (Perez-Jimenez, 2014), they also allow for an alternative means of representing the underlying biocomplexity. By viewing biocomplexity as the length of regulatory events (rather than as a genome size or number of cells), we can approximate how biological processes themselves evolve while still remaining active.
The question arises as to how cybernetic units become ordered in a heterogeneous system. This is similar to the topic of disorder-to-order transition (von Forester, 2003) in the realm of complexity theory. Emergent order through interactions can be understood by breaking down the system into individual behaviors ordered in time. In this formulation, the functional organization cybernetic elements will be ordinal with respect to time, which then imposes order through constraint upon the entire system. This can be thought of as a process of addition, where complexity is added to the emerging system as a series of smaller subsystems. In an embryo specifically, this could is instantiated as a series of timed cell division events. One way to model this type of spontaneous behavior is to look at first-mover dynamics, which is derived from gametheoretic Stackleberg competition (Leitmann, 1978). The computational mechanism we will propose for the emergence of large-scale sets of cybernetic elements, which comprise complex, multicellular embryos, is the Stackleberg leaderfollower game (Hu and Fukushima, 2015). Leader-follower games are generally pairwise, and can be composed of multiple leaders and multiple followers (Sinha et al., 2014). This can occur in as a series of subgames occurring in multiple time-dependent stages, so that each leader-follower subgame can build on the next. A simple demonstration of first-mover dynamics is how the twoplayer competitive game of tic-tac-toe proceeds according to Stackleberg leader-follower rules, and can be mapped to the early rounds of cell division in embryogenesis (Figure 4). While each move in the game (playing a character, dividing a cell) represents a local intentionality, the results (pattern in 3x3 grid, whole embryo) at the global scale are quite unintentional.
In Figure 4A, we can see that X chooses first, and the play is to occupy the center square. Assuming each player makes an optimal move in the context of the previous player's move, the result of this first move is to constrain O into a defensive move. This meta-strategy proceeds until the game ends in a tie, and neither party obtains three characters in a row. First-mover dynamics can also be observed for a dividing embryo in Figure 4B, where division at the one-cell stage produces two cells, which constrain the position and orientation of subsequent divisions. As with the tic-tac-toe example, we make several assumptions about optimality with respect to the size/orientation of cells at each round of division. Figure 4. A diagram of how first-mover dynamics proceed in a game of tic-tac-toe. A: a step-bystep analysis of the first-mover tic-tac-toe-game. Gray squares are the number of possible moves in the next turn. B: a step-by-step analysis of how first-mover dynamics proceeds in a dividing embryo, with a cell's position in a lineage tree shown by sequential numbers inside the cell.
Even though most applications of Stackleberg games have involved non-cooperative decision-making (Pita et al., 2010), first-mover dynamics can be achieved through non-cognitive decision-making . In the tic-tac-toe example, the outcome is determined solely by the first mover's initial move, and in an evolutionary context can be explained through maximization of individual fitness and evolutionary constraint. In fact, each move tends to constrain the subsequent move, and leads to a generally progressive decrease in the system's degrees of freedom over time. An analysis of the tic-tac-toe game in Figure 4A is shown in Figure  5.
In conclusion, we can consider whether the Vehicles model is suitable for understanding evolutionary variation. In Braitenberg (1984), a variety of Vehicle phenotypes are decomposed using a syllabic network. Syllabic networks are constructed by using common features within a description for each vehicle, and uncover commonalities at the level of syllabic similarity. As each vehicle added to the analysis increases in descriptive complexity, new nodes are required to represent this complexity. Such a result is one of an expanded network, both in the number of connections to a single node and overall connectivity of the network. This is analogous to the production of new cybernetic elements during the course of evolution. An expansion in both the number and variety (Ashby, 1958) of cybernetic elements is necessary to allow for increases in complexity without sacrificing functional losses. Expansion of systemic size and complexity may also provide clues as to the level of selection and fitness criteria that allow for regulatory systems to expand while remaining functionally integrated. Figure 5. A graph showing an analysis of the tic-tac-toe example of first-mover dynamics ( Figure 4A), plotted as the number of possible moves per turn versus the identity of the mover (e.g. X or O) per unit time and sequential turn.

Embryogenesis as a Quasi-Reversible Process
Another component of the cybernetic embryo is to understand embryogenesis as a reversible process. Not only do we want to understand what meta-processes aid in constructing the embryo, we also want to identify the essential nature specific combinations of events that guide an embryo through development. In other words, embryogenesis can be both highly deterministic (reversible) or highly contingent (non-reversible). While the reversibility of embryogenesis and similar processes is limited by thermodynamics (Yuan et al., 2015), it may be reversible enough to discover key features through a process of reverse engineering. While we do not as of yet have a clear methodology, experiments utilizing genetic algorithms (Holland, 1992) can inform our understanding of the complexity and reversibility of embryogenesis.
Genetic algorithms are often used to optimize complex processes in a high dimensional and often rough phase space (Kauffman, 1993) that mimics the combinatorial complexity of biological evolution. These spaces resemble the epigenetic landscapes of development, and thus are also applicable to understanding the relationships between development and evolution (Cortez et al., 2017). It encompasses such biological counterparts as genes, genotypes, and phenotypes that are evaluated by some fitness function. Mutations and genotype mating with gene cross-over are also used to produce fitter offspring within an evolving population of genotypes. Fitter offspring are selected to populate future generations. Genetic algorithms have been found to be capable of finding solutions within large and rough feature spaces.
One striking example of the power of artificial evolution as an optimizer comes from Levin's work on Flatworm regeneration (Lobo et al., 2012;Lobo and Levin, 2015;Graeber, 2016). In this case, an optimization procedure was used to discover a model that would explain regeneration and associated phenotypic characteristics of the flatworm. Finding a series of developmental pathways to describe the cascade of events that leads to a successful recapitulation of the phenotype is an exercise in evaluating whether or not the developmental process is reversible. If so, the optimization process should be relatively easy. This can be done in silico in a way that mimics the mechanisms and processes of biological systems (Ray, 1991;Adami and Brown, 1994). Lobo and Levin (2015) used a genetic algorithm to find the head-trunk-tail configuration in flatworm in order to simulate then length of time it takes to recapitulate the regenerative phenotype using no particular pathway. This required six billion simulated experiments and 26,727 generations of models to converge on the desired result. This requires a very large amount of time in evolutionary terms, so it is clear that many of these pathways must be buffered against by the genotype. This approach also works to characterize connectomics in small nervous systems such as that of the nematode Caenorhabditis elegans (302 neurons). In one case (Portegys, 2015), a hybrid genetic algorithm was used to weight synaptic connections in the C. elegans nervous system. There are over 3,000 synaptic connections between these neurons. Synapses are weighted such that the activation states of source neurons affect target neurons variably. These weights could also be optimized could be optimized to produce arbitrary input-output sequences, suggesting possible weighting schemes for the actual synapses to produce observed behaviors. This degree of complexity could be achieved using metamorphic code, which is self-editing but preserves the original function by adding code to the program that accomplish a sub-function through alternative means (Rad et al., 2012).

Conclusions
The embryo is a complex system that increases its phenotypic complexity during a period of great differentiation and specialization. As a process, embryogenesis seems to be both highly regulated and goal-directed towards entropic optimality. Both of these properties can be described in the language of cybernetics. In cybernetic terms, goal-directedness is defined as systems with coherent states that are maintained in the face of perturbation (Heylighen and Joslyn, 2001). In this chapter, we demonstrate this by addressing two questions posed in the Introduction. The first question is answered in Section 6 in the form of black box mechanisms. The second answer, however, is more elusive. Based on our argument in this paper, the mechanism is selection on the various degrees of freedom (or developmental pathways) from single cell to fully developed organism. Selection itself is imposed by thermodynamic entropy and information, the details of which are the subject of future work. Nevertheless, a formal structural model of selection has now been proposed as a set of regulatory mechanisms resulting in patterns such as feedback and crossinhibition (for applications to neuronal systems, see Rutishauser et al., 2015).
Thus, patterns of regulation observed in embryogenesis can be said to be a product of natural selection. Natural selection for goal-directed behavior has a number of consequences. Selection for stability in developmental cognitive systems suggests that positive feedback-driven phase transition serves as the predominant mechanism for the emergence of goal-directedness and even agency (Kelso, 2016). While this example relies on the existence of a formalized nervous system, there is increasing evidence that from a purely information-processing standpoint, such computations can be carried out by other biological circuits and systems . With respect to embryogenesis, two processes in particular are highly similar to formalized neural systems: bioelectric processing between cells during morphogenesis (Levin et al., 2006) and mechanisms that guide cellular behavior on the basis of past experience (Wang et al., 2009;Balaskas et al., 2012). In other multicellular systems such as the slime mold Physarum polycephalum, simple decision-making can occur through cooperativity, positive feedback, and the employment of behavioral heuristics at the organismal level (Reid et al., 2016).
Another way of viewing evolutionary transitions specific to the complexity of embryos involves the concept of metasystem transitions (Heylighen et al., 1995). Metasystem transitions can be defined as emergent processes in associative systems, and consist of composite elements that undergo sudden organizational changes (Francois, 1995). More specifically, metasystem transitions are related to the creation and reorganization of both hierarchical and heterarchical subsystems (Last, 2015). The expected signature of metasystem transitions resembles long periods featuring stable hierarchies and heterarchies interrupted by rapid reorganizations of the relationships between cells (Wake, 1983). In the context of embryos, the "metasystem" allows for collective populations of cells to reorganize both during development and over evolutionary time.
One of our arguments made throughout this chapter involves evolutionary and thermodynamic factors working in concert to direct adaptive processes related to function. This has a conceptual relationship to cybernetics that might inspire future work. Watanabe (1966) has made the connection between Hamiltonian landscapes, the law of requisite variety, and the cybernetic regulation of emergent complexity. Specifically, he cites the inverse H-law as a potential mechanism for the regulation of complexity in emerging multicellular systems. The inverse H-law allows for anti-entropic behavior to guide nonequilibrium systems, or systems that have been perturbed away from equilibrium. This is similar to the argument of Wissner-Gross and Freer (2013), who argue that anti-entropic system behavior is not only expected of the process of entropy maximization in nonequilibrium systems, but leads to the emergence of life and increases in intelligence.
In the case of the cybernetic embryo, the key is how both the regulatory landscape and the epigenetic landscape deal with variation. In the case of the regulatory landscape, the discrete dynamics of a system can account for significant variety if system states are ergodic with respect to a system's dynamics. In other words, control is more or less assured as long as the variation can be anticipated. This might occur if a regualtory system adds components as they are needed in a stepwise fashion (for an example, see Rebeiz et al., 2009). The epigenetic landscape poses a different problem, and requires an additional buffering mechanism (Pujadas and Feinberg, 2012) to "match" the variety encountered as a result of environmental perturbation in a manner consistent with the cybernetic matching and proportional laws of Kompala et al. (1986). In this case, regulatory elements that are in normal times redundant must be evolved to enable ergodicity when variety is encountered. It is likely that if the regulatory system is question encounters unanticipated perturbations, then the addition of regulatory elements must be redundant (addition of elements in a non-stepwise fashion).
Aside from the role of evolutionary and thermodynamic processes in maintaining regulation, there is also a case to be made for the emerging embryo as an active modeler of both evolutionary processes and its own dynamics. In Ashby (1963), the notion of a system modeling its own behavior through approximation is relevant, and can be updated to reflect what we see in embryogenesis. Returning to the notion of multicellular systems emerging from single-celled systems, we see that embryos are governed by both spatial scale and the flow of time. In the case of spatial scale, we have so far discussed the bottom-up aspects of biological complexity. However, there is also an element of "wholeness" (Bohm, 1980) in the embryo, particularly when one considers the embryo as a single, discrete entity (Gordon and Gordon, 2016). The flow of time is also a particularly interesting problem that goes beyond the conventional aspects of dynamic regulation. In terms of regulation, one might view the embryo and its constituent parts as having a retrocausal temporal structure (Krasnikov, 1997). Specifically, the component parts (genes, proteins, and individual cells) might change in response to selective pressures at the level of the whole embryo. This rather specialized view of multilevel selection (Krupp, 2016) would allow for changes to unfold in the development and evolution of future embryos, but also allows for selection to act upon lower organizational levels and their own regulatory requirements.
In the case of embryogenesis, the context of developmental contingency could determine whether the process of embryogenesis is either hierarchical or heterarchical. Hierarchical structures (such as found in the lineage tree of the nematode C. elegans) may emerge as a response to high degrees of environmental perturbation. This would ensure consistency in the face of unexpected variation. By contrast, heterarchical structures (such as found in Vertebrates) might emerge in systems that have additional regulatory mechanisms to buffer unexpected variation. The cybernetic embryo also provides a critical link between emergence, evolution, and adaptive biological processes. For example, there are also connections that link interactivity and regulatory motifs that are worthy of further study, such as how tic-tac-toe played in Stackleberg equilibrium resembles a futile cycle. In this case, the set of resulting interactions produced thermodynamic entropy without a clearly advantageous outcome for each competitor. While the first mover in a pairwise interaction constrains the second mover, playing the game optimally results in a cancellation of initial advantage that resembles the lack of production in a futile cycle. However, these interactions do result in one useful outcome: structure potentially used as a scaffold for increasing complexity.
There are many potential criticisms of our approach, particularly as it simultaneously covers much intellectual ground while also claiming great generalizability. The abstract nature of our approach, which is partially reliant on biosemiotic and other symbolic approaches, is also quite unconventional. There have been systems biological critiques of semiotic approaches, particularly as to whether or not natural phenomena can be reduced to natural (or artificial) linguistic forms (Emmeche and Hoffmeyer, 1991). However, there are also two advantages to using such an approach. The first is that a general structuralist approach is a useful heuristic in dealing with the holistic nature of phenotypes (Emmeche, 2011). Additionally, emergentist approaches also help us understand emergent processes (Emmeche and Hoffmeyer, 1991) and act to countervail genecentric approaches (Alicea and Gordon, 2014). Along with these advantages, our approach provides a formal set of links between the evolution of development and regulatory mechanisms provided by the cybernetic embryo framework is a crucial component to reverse engineering the embryo defined as a multicausal, cascading, nested decision-making system (Gordon and Stone, 2017). In this context, differentiation is the result of the arrival of a networked decision. Feedback signals, across space and time, are a metaphorical 'spooky glue' holding the overall physiochemical system of evolution and development together. Future work might focus on building meta-models of this process that include multiple scales of biological organization.


s2fieldsofstudy: Biology Medicine

text:
Somatosensory imagery induces topographically specific activation patterns instrumental to fMRI-based Brain Computer Interfacing

For patients with ‘locked-in’ syndrome (LIS), brain-computer interfaces (BCIs) based on hemodynamics are a promising tool to regain communication. Visual and motor imagery are oftenused information-encoding strategies despite frequent visuo-motor impairments in LIS. However, evoking imagery not grounded in recent sensory experience is very challenging. Therefore, somatosensory imagery might be more suitable, as somatosensation is preserved even in late-stage LIS. Using 7T fMRI in healthy subjects, we find that somatosensory imagery is associated with topographic decoding-weight patterns in primary somatosensory cortex and show high-accuracy decoding with little pre-training. We immediately implement these novel findings to drive brain-based communication, providing a promising proof-of-concept for a novel type of communication BCI with high clinical potential. Our new BCI-control strategy paves the way for several future applications. Next to brain-based communication, it might be employed for neurofeedback-based training to restore somatotopy in stroke and chronic pain patients suffering somatosensory deficits.

Introduction
suitable for online BCI-based communication; (2) explore to what extent its classification is driven by somatotopic information. For this proof-of-concept study, we exploited the high signal to noise ratio of ultra-high field 7T functional magnetic resonance imaging.

Results
Ten healthy participants were instructed by auditory cues to perform alternating somatosensory imagery of the hand or foot (figure 1) while ultra-high field 7T fMRI data were collected. Each participant completed five (participant 7) to six somatosensory-imagery classifier-training runs. Seven participants also completed one or more answer-encoding runs in which they answered a binary question using imagined somatosensory stimulation to the right hand to encode "yes" and imagined somatosensory stimulation to the left foot to encode a "no" answer.

Suitability of somatosensory imagery for BCI
After establishing that all participants subjectively experienced somatosensory imagery (see methods section 4.2), we went on to test the suitability of this strategy to induce activation that can be reliably decoded for the purpose for BCI communication. We hypothesized that, using a machinelearning algorithm, the contents of somatosensory hand or foot imagery can be accurately decoded on the single-trial level with similar accuracy for each imagery class, with accuracy increasing when increasing numbers of runs are used to train the classifier (hypothesis 1a). Our second hypothesis (1b) was that somatosensory imagery would yield accurate decoding of answer runs in simulated real-time as well as during a proof-of-concept online brain-based communication session.
Testing the first hypothesis (1a), we found that significant above-chance classification accuracy (p<0.05, based on permutation testing) was achieved in eight out of ten subjects with one training run, and for all subjects for two or more trainings runs (figure 2). The average proportion correct classifications increased with the number of runs included for training (figure 3), up to 0.82 (sd= 0.12) when using five runs (four in one subject) for training and one run for testing. The effect of imagery class and the amount of training runs was tested in an imagery class (2, foot or hand) by training runs To test the second hypothesis (1b) and provide a proof-of-concept for brain-based communication, we used the answer runs that were available from seven subjects (17 runs in total, 7 'no', 10 'yes'; table 1). A simulated real-time analysis revealed no significant effect of the number of runs used for training (F(5, 30) = 1.470, p = 0.229, figure 4). Using all training runs, the average proportion correct for these seven subjects was 0.86 (sd = 0.17), mounting to 0.92 (sd =0.08) when disregarding one subject who, after the session, indicated he misunderstood the instructions for the answer runs. In the proof-of-concept online communication session with the last subject all answers to the questions were accurately decoded.

Assessment of somatotopy
Secondly, we wanted to assess whether classification success would be based on somatotopic information. We hypothesized (2a) that classification accuracy would be higher in primary somatosensory cortex, where somatotopy is most pronounced, compared to secondary somatosensory cortex and (2b) that the discriminative weights obtained from the support vector machine (SVM) would show a somatotopic distribution.
Next, we tested somatotopy of the discriminative SVM weights (hypothesis 2b). Visual inspection of discriminative SVM weight maps revealed a pattern that roughly corresponded to the expected somatotopy in four out of ten subjects (2, 8, 9 and 10; figure 6). A weaker somatotopic pattern was observed in one subject (4). A reversed pattern was noticed in two subjects (3 and 6), and a mixed pattern was seen in three final subjects (1, 5 and 7). In a more quantitative evaluation of the extent to which there was a somatotopic distribution in the discriminative SVM weights, we tested whether there was a distinction in the median coordinates of the voxels associated with the 200 largest positive (foot) weights and the largest negative (hand) weights. We focused on the x-coordinates only, because somatotopy is expected to be most pronounced on the medial-lateral axis of the brain. A Wilcoxon signed rank test revealed that there was a significant within-subject difference (z = 2.09, p = 0.0365) between the median x-coordinates of the voxels associated with the 200 largest positive (foot) weights and the median x-coordinate of the voxels associated with the 200 largest negative (hand) weights.
The median x-coordinate across subjects was 8.50 (sd = 14.5) for the 200 largest positive weights and -21.5 (sd = 18.4) for the 200 largest negative weights (figure 7). Also, we found a significant correlation between the difference in median x coordinates and the overall classification accuracy, i.e. subjects with a larger value for the foot-hand median x-coordinate difference (larger 'hand-foot distance') also showed a higher classification accuracy (Spearmann's rho .948, p = 0.000).

Discussion
The current study (1) investigated whether somatosensory imagery is a suitable encoding strategy for a hemodynamic (fMRI-based) BCI, aiming to provide a proof-of-concept of its suitability for online BCI-based communication and (2) tested to what extent its decoding is based on somatotopic information. In accordance with our first set of hypotheses, we found that the contents of somatosensory imagery could be decoded above chance in each of our participants with similar accuracy for hand or foot, and that decoding accuracy increased with an increasing number of training runs. We also found that answers to binary questions encoded by seven healthy participants using the somatosensory imagery BCI could be successfully decoded in simulated real-time as well as during a proof-of-concept online communication session tested in one subject. In accordance with our second set of hypotheses, we found that decoding accuracy was significantly higher in primary somatosensory cortex compared to secondary somatosensory cortex, and that discriminative weight patterns generated by healthy participants through somatosensory mental imagery showed a somatotopic distribution.

Suitability for use in fMRI Brain-Computer-interfacing
The average classification success was 82% for the offline leave-one-run out cross-validation procedure on the training runs, and 92% for the classification of the answer runs in simulated realtime (86% when including a subject who indicated afterwards that he misunderstood the instructions).
The fact that the classification success was higher for simulated real-time classification of the answer runs could be explained by the additional training run that was used for training the classifier in simulated real-time, whereas in the offline leave-one-run out procedure one run necessarily had to be left out for testing. However, the evaluation of the effect of additional training runs in the offline condition suggested that the beneficial effect of adding training runs levelled off after four runs.
Another possible explanation might lie in the slight differences in the analysis pipelines for the simulated real-time and offline analyses regarding slice-scan-time-correction and data transformation (Talairach transformation and downsampling to 2mm 3 offline versus smoothing with 4mm FWHM Gaussian kernel in native space in simulated-real-time). On the other hand, the increased accuracy for the answer run classification might also be related to enhanced imagery quality during the answer run, due to extended practice, focus on a single imagery strategy (in contrast to the alternating imagery strategies in the training runs) and possibly increased motivation because of the social interaction component of answering a genuine question.
Taken together, these findings inspire confidence that somatosensory imagery will be a viable additional control strategy for BCI communication procedures. Somatosensory imagery might be especially suited for locked-in patients with reduced visual or auditory faculties who lack muscle control. These patients might not be able to evoke stable imagery involving their most affected modalities, as they lack recent memories of sensory experiences which can increase the level of detail and vividness of mental representations generated by imagination (Hassabis & Maguire, 2009). The current study was performed at 7T, but given the relatively large distance between hand and foot representations in primary somatosensory cortex similar results could be expected at a field strength of 3T which is more commonly available in clinical settings. Our novel somatosensory imagery control strategy might also be applied in binary fNIRS-based communication BCIs. Of course, taking into account the considerably lower spatial resolution of fNIRS, optode placement and data analysis procedures would have to be optimized.

Topographic nature of somatosensory imagery
Our results indicate that decoding accuracy was higher within the primary than in the secondary somatosensory cortex mask, and that it is possible to induce somatotopic discriminative weight patterns by mental imagery. The median x-coordinates of the subjects with highest offline classification accuracies roughly correspond to values reported in the literature for finger (Kolasinski et al., 2016;Nelson & Chen, 2008;Pfannmoller, Greiner, Balasubramanian, & Lotze, 2016;Sanchez-Panchuelo, Francis, Bowtell, & Schluppeck, 2010;Schweizer, Voit, & Frahm, 2008) and foot representations (Akselrod et al., 2017; figure 7). The slightly more medial imagery-related weight coordinates could be due to the fact that, in the current study, imagery involved the whole hand whereas stimulation is usually delivered to the fingers only. For the foot, imagery coordinates were more inferior than reported stimulation coordinates. Similar observations of corresponding sites for imagery and stimulation in primary sensory cortices have also been made in previous studies in the somatosensory domain (de Borst & de Gelder, 2017;Newman et al., 2005;Schmidt et al., 2014;Wise et al., 2016;Yoo et al., 2003), however, to our knowledge, the current study is the first to report a somatotopic pattern of discriminative weights in the (primary) somatosensory cortex during imagery.
Previous studies have reported topographic recruitment of visual cortex during imagery (Klein et al., 2004;Slotnick, Thompson, & Kosslyn, 2005;Thirion et al., 2006) The presence of a somatotopic pattern in the discriminative weights, that is a more medial and right sided location for the highest positive weights corresponding to the imagery class 'foot' and a more lateral and left sided location for the highest negative weights associated with 'hand', was associated with higher classification success. This suggests that subjects might increase classification success by increasing the level of somatotopy in their imagery-induced brain activation patterns. If this finding can be extrapolated to subjects with maladaptive somatosensory representations, we speculate that the current BCI might be suitable to use for training the somatotopy of these representations. Somatosensory dysfunction and aberrant plasticity have been linked to pain severity in chronic pain (Di Pietro et al., 2013;Di Pietro, Stanton, Moseley, Lotze, & McAuley, 2015;Flor, Braun, Elbert, & Birbaumer, 1997;Kim, Kim, & Nabekura, 2017;Pleger et al., 2005;Wrigley et al., 2009), but see Makin et al. (2013) for an alternative view). Neurofeedback based on somatosensory mental imagery could be especially suitable for clinical interventions in patients for whom actual somatosensory stimulation is difficult or impossible due to their physical limitations or technical limitations of the scanner environment.
Somatotopic patterns could in theory also be evoked using motor imagery or even actual movement. However, our subjects were carefully instructed and trained in a mock scanner to exclusively use somatosensory imagery, and it was emphasized that they should not move during imagery runs. Subjects' compliance to these instructions is supported by the fact that we did not observe hand or foot movement during training or scanning. Also, during the training runs, subjects did not receive any feedback or reward based on their performance, so there was no incentive to try out non-somatosensory imagery strategies. In line with this, subjects' reports indeed indicated that they used purely somatosensory strategies (supplementary table 1).
Based on the current data, we cannot establish a link between classification accuracy and the use of a specific somatosensory imagery strategy: The two subjects with largest classification success used touch and warmth imagery, as did the two subjects with the lowest classification success. It would be interesting to investigate whether imagery targeting different somatosensory submodalities induces different patterns of activation in primary somatosensory cortex. However, such pattern differences might be most apparent within a limb-specific subregion, whereas the general somatotopy of hand and foot is most likely preserved across somatosensory submodalities (Sur, Wall, & Kaas, 1984).
We instructed our subjects to stick with the same strategy for hand and foot. Whereas some subjects only reported using imagery based on one somatosensory submodality (e.g. vibrotactile stimulation), others used imagery of more naturalistic somatosensory stimulation involving a combination of somatosensory submodalities. Such imagery, e.g., combining warmth and touch might increase vividness and evoke stronger activation. Work by Kosslyn and Thompson in the visual domain suggests that the level of vividness and detail is related to primary visual cortex activation during imagery (Kosslyn & Thompson, 2003). We therefore propose to leave the choice of strategy to the participant, and have them choose the strategy that subjectively evokes most vivid and detailed imagery and that they feel most comfortable with.
There are some limitations to this study. Somatosensory imagery involved only two contralateral body parts, which are relatively far apart in physical and somatotopic space, and whose representations reside on contralateral sides of the body. Thus, classification could have been based solely on less specific imagery related to body side. However, the fact that classification success was highest within the primary somatosensory cortex mask and linked to the degree of somatotopy in the discriminative weight maps indicates that subjects used a limb-specific somatosensory imagery strategy rather than a body-side related strategy. This is also in accordance with the introspective description of imagery strategies reported by the subjects. However, note that the question whether the brain signal is body-part or body-side specific is irrelevant in the context of the suggested twochoice communication BCI application, as long as it can be used at will.
A final limitation of this study is that only healthy subjects without neurological impairments participated in the current study. Patients might have shorter attentional spans leading to more 1 0 variable imagery quality. However, we found higher classification rates during answer runs than during training runs, possibly due to increased motivation. This effect might positively extrapolate to patients in a BCI communication setting.

Conclusion
Taken together, we conclude that somatosensory imagery is a successful informationencoding strategy for hemodynamic BCI-based communication with a high clinical potential. The current results pave the way for future applications, such as brain-based communication for patients with locked-in-syndrome, and neurofeedback-based training to support the rehabilitation of somatosensory representations in stroke and chronic pain patients suffering from somatosensory deficits.

Participants
Ten neurologically healthy participants (six males, mean age = 29 years, sd = 5 years, one left-handed) were recruited among students and staff members of the Faculty of Psychology and

Procedure and experimental design
Two to five days before the actual scanning took place, subjects participated in a training session in which they performed one or more imagery runs in a mock scanner. The purpose of this training session was to instruct them about the task and possible imagery strategies, and to encourage them to try several imagery strategies and select the one that worked best for them, evoking the most stable and vivid somatosensory experience, as assessed by introspection. Two imagery strategies were suggested as possible starting points. The first strategy was to imagine someone gently touching the inner surface of the hand or foot. The second strategy was to apply techniques from autogenic training, imagining a heavy and warm feeling in the hand or foot. The session lasted at least 30min or until the participant indicated (s)he had found a subjectively effective imagery strategy.
During the scanning session, the participants' heads were stabilized with foam padding to minimize head movement. They were instructed not to move and to keep their eyes closed during the 1 functional runs. They received auditory instructions indicating when they had to start imagining touch to the right hand ("hand") and the left foot ("foot"), and when they could stop ("rest"). Each participant completed five (participant 7) to six somatosensory-imagery classifier-training runs. Every classifier-training run contained a pseudorandom sequence of nine left-foot (LF) and nine right-hand (RH) somatosensory mental imagery trials of 18s duration, each followed by a 16, 18 or 20s rest period (figure 1). Seven participants also completed one or more answer-encoding runs in which they answered a binary question using imagined touch to the right hand to encode "yes" and imagined touch to the left foot to encode a "no" answer. Each answer-encoding run included five tactileimagery trials (encoding the same answer). Timing of the imagery trials and rest intervals was identical to the training runs. If possible, each question was followed by two answer runs. In the first answer run, subjects were asked to encode the actual answer to the question using the appropriate imagery strategy. In the second answer run, they were asked to encode the opposite answer (in order to obtain approximately the same number of "yes" and "no" answers-encoding data). The correct answers to the questions were unknown to the researcher responsible for data analysis until after completion of the offline analyses. In one subject, the answer runs were also analysed online and the results were reported back to the participant right after the run.
After scanning, participants were asked to fill out a questionnaire in which they had to rate the experienced quality of their somatosensory imagery (i.e., whether it was a robust mental experience) and the extent to which it evoked a tactile sensation (i.e., to what extent did they experience it as inducing tactile input coming from the target limb) on a ten-point scale (one being the worst/least applicable, ten being the best/most applicable). They also rated the percentage of time that they were successful at maintaining the imagery and described the particular imagery strategy they used.
All subjects reported successful somatosensory imagery. The type of strategies used ranged from imagining a warm feeling (four subjects), imagining vibratory stimulation to the finger or toe (two subjects), to imagining touch (massage, brush strokes, somatosensory sensation; four subjects; supplementary table 1). On a scale of 1-10, the average score subjects gave for the quality of their hand imagery was 6.9 (sd = 0.65) and of their foot imagery 6.6 (sd = 0.49). The average score in response to the question to what extent they experienced a tactile sensation during imagery was 5.3 (sd = 2.3) for the hand and 4.9 (sd = 1.9) for the foot. The reported average percentage of time that somatosensory imagery worked during the hand imagery epochs was 75% (sd = 4.5%, N=6) and 73% (sd = 8.2%, N=6) for the foot imagery epochs. Paired t-tests did not reveal any significant differences between foot-and hand-related ratings at p<0.05.

MRI data acquisition
Data were acquired at Maastricht Brain Imaging Centre (Maastricht University) on a 7T scanner (Siemens Healthcare, Germany) with a 32 receiver channel (Nova) head coil.
The scanning session included the acquisition of an anatomical scan (8min) and five to six classifier training runs (12min/run). In seven subjects, there was sufficient time to end the session with one to seven additional answer-encoding runs (1 answer run: 2 subjects, 2 answer runs: 4 subjects, 7 answer runs: 1 subject) in which the participant answered one or more binary questions using somatosensory imagery (4min/run).
In the functional runs, a gradient echo sequence was used with 78 slices of 1.5mm thickness, 1.5mm in-plane resolution, 198 mm x 198 mm field-of view, TR/TE = 2000/26.1ms, 60 degree flip angle, multiband acceleration factor 2 and GRAPPA acceleration factor 3 (340 volumes for the training runs, 96 volumes for the answer runs). In each session, five volumes were acquired in opposite phase encoding direction to be used for (offline) correction of EPI distortions.

Data analysis Preprocessing for offline multivariate analyses
Offline data analysis was performed using Brainvoyager QX 2.8.4 (BrainInnovation, Maastricht, the Netherlands). The anatomical scan was down-sampled to 1mm 3 . The first two volumes of each functional run were discarded. The remaining volumes were adjusted based on slicescan time, realigned to the first functional image of the session using rigid body transformations, temporally filtered to remove linear trends and nonlinear periodic signals with a frequency below 0.01Hz, corrected for EPI distortions, coregistered to the anatomical scan interpolating to 2mm 3 isotropic, and normalized to Talairach space. This procedure includes the definition of subject specific anatomical landmarks (AC [anterior commissure], PC [posterior commissure] and the borders of the cerebrum) that are used to rotate each brain in the AC-PC plane followed by piecewise, linear transformations to fit each brain in the common Talairach 'proportional grid' system (Talairach & Tournoux, 1988).
The multivariate analyses were limited to the voxels included in anatomical masks. Masks were obtained by creating anatomical regions of interest (ROIs) in Talairach space based on the Jülich probability maps for the primary (PSC3a, PSC3b, PSC2 and PSC2;Geyer, Schleicher, & Zilles, 1999;Geyer, Schormann, Mohlberg, & Zilles, 2000;Grefkes, Geyer, Schormann, Roland, & Zilles, 2001) and secondary somatosensory cortex (SIIOP1, SIIOP2,SIIOP3 and SIIOP4; Eickhoff, Amunts, Features (t-values per voxel) were extracted by linearly fitting each voxel within the mask with a design matrix consisting of a constant term (intercept) and a modelled hemodynamic response, within a temporal window of 30s, spanning from one volume (2s) before the auditory cue to 14 volumes (28s) after the auditory cue. The modelled hemodynamic response was obtained by convolving a canonical hemodynamic response function (hrf) with a box-car predictor whose value was one in the volumes corresponding to the imagery blocks, and zero otherwise.

Effect of the amount of training data and imagery class
To test the hypothesis that the amount of training data would increase decoding accuracy independent of the imagery class the number of runs used for SVM training was increased from one to five, each time testing on one of the remaining runs, using all possible combinations. Features were extracted from the largest, most inclusive anatomical masks of the primary and secondary somatosensory cortex. To test single-subject significance, permutation tests were performed. The edge corrected and logit transformed classification accuracies were subsequently entered in a two-way repeated-measures ANOVA in with imagery class (2) and number of training runs (5) as within-subject factors. To test whether accuracy increased with the number of training runs, within-subject planned polynomial contrasts were used.

Assessment in the context of BCI communication
To test the suitability of the current approach for (online) communication, a simulated realtime classification was performed using Turbo-BrainVoyager (version 3.2, BrainInnovation B.V., Maastricht, the Netherlands) to emulate the procedure that would be performed during the scanning session. In the last participant, a proof-of-concept online communication experiment was performed in which the answer runs were also decoded online during the scanning session using an SVM based on data from all available training runs.
The online and simulated real-time analyses were based on slice data in native resolution (1.5mm 3 ), which were 3D-motion corrected using the standard TBV incremental procedure spatially realigning each volume to the first recorded volume by rigid body transformation. The lowest probability (largest) S1 and S2 mask was back-transformed from Talairach space into native space.
After spatially smoothing the functional data with a 4 mm full -width -at -half -maximum (FWHM) kernel, t-values were computed by fitting a general linear model with the same predictors as used in the offline procedure, within the same temporal window of 30s, spanning from one volume (2s) before the auditory cue to 14 volumes (28s) after the auditory cue. T-values from the S1 and S2 mask were used as features for SVM training and testing using LIBSVM (2000-2009 Chih-Chung Chang andChih-Jen Lin).
In the simulated real-time analyses, the effect of the number of training runs was tested by using a classifier based on an increasing number of runs, starting from the run closest to the answer run to emulate the temporal adjacency that is characteristic for the online situation, each time adding a preceding run for each additional level for SVM training, and subsequently testing on the available answer run(s). The resulting edge-corrected and logit transformed accuracies were entered into a repeated-measures ANOVA with number of training runs (1-6) as within factor was used to statistically assess the effect of the amount of training data.

Assessment of somatotopy Comparing accuracy in primary and secondary somatosensory cortex masks
The effect on decoding accuracy of the anatomical region chosen for feature extraction was evaluated by repeating the cross-validation procedure with five runs (four for subject 7) for training and one run for testing with features extracted from the anatomical masks of the primary (PSC3a, PSC3b, PSC2 and PSC2) and secondary (SIIOP1, SIIOP2, SIIOP3 and SIIOP4) somatosensory cortex at different probability levels (Eickhoff, Amunts, et al., 2006;Eickhoff, Schleicher, et al., 2006;Geyer et al., 1999;Geyer et al., 2000;Grefkes et al., 2001). The probability level (a percentage) of a mask for a particular anatomical area indicates that all voxels in that mask were classified as part of that anatomical area in at least the corresponding percentage of subjects in the database used to construct the anatomical probability map. For example, an anatomical mask with a probability level of 30% contains voxels that were classified as belonging to that particular anatomical region in at least 30% of the subjects in the database used to construct these anatomical probability maps. The higher the probability level, the smaller the mask size (supplementary table 2). To test single-subject significance, permutation tests were performed.
Accuracy values were entered in a two-way repeated-measures ANOVA, with anatomical region (2, S1 or S2) and probability level (8, 10-80%) as within-subject factors, to test whether there was a difference between S1 and S2 in terms of decoding accuracy, and whether this was possibly affected by the probability level, which was further assessed using within-subject planned polynomial contrasts. For the ANOVA, only probability levels of 80% or lower were taken into account which corresponded to masks with more than 500 voxels, to ensure stable classification results.

Assessment of somatotopy in discriminative SVM weights
The extent to which the classification was based on somatotopically different weight patterns in primary somatosensory cortex was assessed using the average discriminative weights across the different training splits. Discriminative weights were extracted from single-subject maps (obtained when testing the overall single-subject classification accuracies) that represented the average weights across the different training splits, using SVM classifiers trained on features from the most inclusive anatomical mask, with all possible combinations of 4/5 runs for training. Because the representations of hand and foot are most clearly separated on x-axis (left-to right), we computed the subject-specific median x-coordinates for the 200 largest positive (foot imagery) and 200 largest negative weights (hand imagery), and compared them using a Wilcoxon signed rank test. Finally, the Spearman rank correlation was computed between the subject-specific differences in median x-coordinates for foot and hand and the leave-one-run-out classification accuracy revealed using phase-locked analysis and verified using multivoxel pattern analysis and functional connectivity. J Neurosci, 35 (7) Tables   Table 1: Data from simulated real-time analysis of answer runs. No answer-encoding data were obtained for subjects 4, 7 and 9 due to running short of MR scanning time.     Figure 5: Offline single-trial classification accuracies for S1 versus S2. Accuracies were obtained using a leave-one-run out cross-validation scheme based on spatial features extracted from separate masks of S1 (red line) and S2 (purple line) at different probability levels ranging from 10% to 80% . Higher probability levels were excluded because the masks included less than 500 voxels for the S2 mask. The shaded areas indicate the standard deviation for each case.