To swim or not to swim: A population-level model of Xenopus tadpole decision making and locomotor behaviour.

We present a detailed computational model of interacting neuronal populations that mimic the hatchling Xenopus tadpole nervous system. The model includes four sensory pathways, integrators of sensory information, and a central pattern generator (CPG) network. Sensory pathways of different modalities receive inputs from an “environment”; these inputs are then processed and integrated to select the most appropriate locomotor action. The CPG populations execute the selected action, generating output in motor neuron populations. Thus, the model describes a detailed and biologically plausible chain of information processing from external signals to sensors, sensory pathways, integration and decision-making, action selection and execution and finally, generation of appropriate motor activity and behaviour. We show how the model produces appropriate behaviours in response to a selected scenario, which consists of a sequence of “environmental” signals. These behaviours might be relatively complex due to noisy sensory pathways and the possibility of spontaneous actions.


2
The hatchling Xenopus tadpole provides a good place to study decision making and behaviour because at this stage of development (approximately 2 days old, stage 37/38) the nervous system is relatively small (several thousand neurons), the behaviour is simple, and many biological details are known from experimental work.
In this paper, we present a new computational model of the tadpole nervous system that is informed by experimental data and is able to accurately reproduce tadpoles' behaviour in response to input from multiple sensory modalities. We implement the integration of noisy sensory signals in a simple model of interconnected neuronal populations. This model can describe the behavioural switching observed in hatchling Xenopus tadpoles (Roberts et al., 2010). The aim of the model is to clarify the key universal neurobiological mechanisms and theoretical principles that underlie the decision-making process, as well as provide new ideas and hypotheses for experimental testing.
Many prevailing theories regarding decision making postulate that evidence from different sensory modalities is integrated until a threshold is reached and one of several actions is selected (Kristan, 2008;Gold & Shadlen, 2007;Marshal et al., 2012). This integration process is important due to the noisy nature of sensory signals. Our computational model describes the dynamics of behavioural responses to input signals from the "environment". The modelling is based on multiple anatomical and neurophysiological findings regarding initiation of locomotion by trunk skin (Roberts et al., 2014) and head skin (Buhl et al., 2015) stimulation, as well as several other sensory inputs that are known to control tadpole locomotion (Roberts et al., 2010). The model includes two parts: (1) four sensory pathways (touch trunk skin, touch head, light dimming, and press head) and (2) the central pattern generator (CPG) neurons for execution of locomotor actions.
All sensory pathways in the model are organised in a similar way. A neuronal population corresponding to a particular sensory modality processes information from (non-modelled) sensory cells and delivers the result of this processing to a central integrator population. This integrator population is where decision-making and action selection occurs. The CPG populations work under the control of inputs from the integrated sensory signals to generate neuronal activity corresponding to one action from a repertoire including swimming, struggling, and accelerated swimming (Roberts et al., 2010;Li et al., 2015). The output of the CPG is motor neuron spiking, with different patterns of spikes corresponding to different actions.
Many experimental facts are known about hatchling Xenopus tadpoles regarding the different neuronal types that are present, including their anatomy, electrophysiology and synaptic connections. However, the amount of experimental data available is still not enough to produce a detailed single neuron level model of all sensory modalities and behaviours. Therefore, our approach is based on a mean-field (mesoscopic) model of neuronal activity.
The model is formulated as a system of 26 ordinary differential equations, which describe the average activity level in various interacting neuronal populations. The model is symmetrical, with 13 populations on the left-hand side of the body and 13 populations on the right-hand side. We describe the dynamics of population activity using the Wilson-Cowan model (Wilson & Cowan, 1972;Borisyuk & Kirillov, 1992). For each population we describe the dynamics of the average neuronal activity in the population of neurons in response to incoming synaptic inputs from the other populations. We use bifurcation analysis to determine the region in parameter space that corresponds to physiological activity such that the model produces the correct output. For example, for a pair of coupled populations, a region of parameter space where regular oscillations exist can be determined (Borisyuk & Kirillov, 1992). Bifurcation analysis also reveals how parameters can be changed to control the dynamics and switch from one dynamical mode to another.
Simulations show a good agreement between modelling results and experimental measures. For modelling swimming we use a bi-stable regime that exists in the system: a short-term input corresponding to stimulation of the trunk skin on one body side moves the system from a stable equilibrium to generating anti-phase oscillations (alternating left and right activity). This swimming mode exists for some time (if there are no perturbations from the sensory pathway) before spontaneously stopping and returning the system to the resting equilibrium. We show that the head touch sensory pathway can initiate swimming in a similar way.
According to experimental data, struggling behaviour is slower, stronger series of rhythmic trunk flexions seen while a tadpole is grasped by a predator (Roberts et al., 2010) and the corresponding spiking activity is of bursting type. To model struggling we use a dynamical regime of anti-phase envelope oscillations where the slow frequency relates to slow body movement and the fast frequency reflects bursting spiking. We show that this struggling mode in the model appears and exists during prolonged input from the skin touch sensory pathway.
Tadpoles have a pineal eye that is able to sense light. When the light sensed by this eye is dimmed, the swimming frequency increases (Jamieson & Roberts, 2000) which has the effect of causing the tadpole to swim upwards. The model mimics these experimental findings, increasing the activity rate of swimming transiently in response to input from the light dimming sensory pathway.
Experiments with the cement gland sensory pathway reveal that swimming stops when the cement gland is pressed or its mucus pulled (Roberts et al., 2010). A tadpole also stops swimming when it bumps into solid objects like vegetation or the side of a dish. Model simulations show that inhibitory input from the "head press" (cement gland) sensory pathway on one body side stops swimming and returns the system to the resting equilibrium.
We demonstrate that for a selected sequence of events (touch skin, dimming light, predator attack, head bump) reflecting a "natural" scenario (Li et al., 2014), the model correctly reproduces real CPG activity. These neuronal activities correspond to the tadpole's behaviour. Thus, the model can generate proper behaviour in response to external events through the integration of sensory inputs and decision-making.

Repertoire of tadpole behaviours
In this section, we provide a short introduction to the set of possible tadpole behaviours (Li et al., 2014), which are reflected in our modelling.
In the resting state a tadpole is usually attached to some object (water surface, wall or bottom of a dish, etc). After about 20s in this state the tadpole will drift to the bottom of the water and stay there for some time (about 60-90s), and after that spontaneously start to swim (Jamieson & Roberts, 2000). Swimming consists of alternating left-right motor neuron activity at 10-25Hz. If a tadpole is not moving (steady state) then swimming can start by one of the following ways: (1) on trunk skin touch; (2) on head skin touch; (3) on light dimming; (4) spontaneously. Swimming can stop: (1) on press head or (2) spontaneously.
Struggling behaviour is slower, stronger series of rhythmic trunk flexions seen while a predator grasps a tadpole. During struggling, active neurons fire bursts of spikes (Li et al., 2007). Each burst is about 100-200ms long and the intra-burst spiking frequency is up to 245Hz. Struggling can start (Soffe, 1991) in response to prolonged stimulation of trunk/head skin on both sides simultaneously, either when the tadpole is at rest or swimming. Struggling lasts during the stimulation period only and after that the tadpole switches to swimming.
When a tadpole is swimming, dimming the light causes the frequency of swimming to increase, which causes the animal to swim upwards (Jamieson & Roberts, 2000). This swimming with increased frequency lasts for the duration of sensory input caused by the light dimming; if the light level returns to normal then swimming switches back to its normal frequency.

Model formulation
This sections begins with an explanation of the neurobiological details of sensory pathways and CPG neuron populations and their functional mechanisms in the Xenopus tadpole. Using these facts we formulate a model of interactive neuronal populations. Each population that has been included to the model corresponds to a set of neurons of some particular type in the real tadpole. The coupling between populations relates to experimentally measured neuronal connectivity. The dynamics of the average level of activity in each population is described by the Wilson-Cowan equation (Wilson & Cowan, 1972). This equation is a simplification of real neurophysiological processes; however, this model reflects the most important details of dynamics at the population level. We derive equations for both sensory pathway and CPG networks and describe how the sensory pathway delivers information from external inputs to CPG populations to produce proper motor neuron activity.

Biological motivation for the model
In this section, we specify some neurobiological details of the tadpole nervous system which are important for the model formulation (Li et al., 2014). We consider two parts of the nervous system: sensory pathways (SP) and the central pattern generator (CPG).

Sensory pathways.
We consider the following four sensory pathways ( Fig. 1): (1) Touch trunk skin (TS). This pathway includes Rohon-Beard (RB) cells, which are distributed along the rostro-caudal body dimension (about 100 neurons per side). In response to touching the trunk skin, several nearby RB neurons will each generate a single spike (Roberts et al., 2014). These spikes are propagated by the sensory pathway: dorso-lateral commissural and ascending neurons (dlc and dla, respectively). We assume that there is a population of interneurons (xINs) which receives signals from both dlc and dla neurons and tINs (see below). The population of xINs integrates incoming signals and delivers the excitation to dIN and dINr populations in the CPG.
(2) Touch head (TH). This pathway includes about 70 trigeminal sensory touch receptors (tSt) on each body side, which generate a single spike on touch. These sensor excite the neurons of the trigeminal nucleus (trigeminal inter-neurons; tINs) (Buhl et al., 2015). We assume that the tIN population excites the same integrating xIN populations as the trunk skin pathway (Fig. 1). The output from xIN populations builds up activity in the descending interneurons (dINs) of the CPG, potentially resulting in a decision to start swimming. (3) Light dimming (LD). This pathway includes the photoreceptors of the pineal eye, which innervate pineal ganglion cells (pgc), and diencephalic/mesencephalic descending (D/Md) neurons on both body sides (Jamieson & Roberts, 1999). We assume that the output of the LD pathways converges onto dINs directly, bypassing the xIN populations. In response to light dimming, a tadpole starts swimming or, if already swimming, accelerates during the period of light dimming (Jamieson & Roberts, 2000). (4) Press head (PH) inhibitory pathway. Cement gland receptors (tSp) produce spikes on when pressure is applied to the head and innervate the inhibitory population of mid-hindbrain reticulospinal (MHR) neurons on both body sides (Perrins et al., 2002;Li et al., 2014). This inhibition converges to dINs and other CPG neurons (including motor neurons) to stop their activity. Figure 1 summarises the connections between populations in the sensory pathways. This connectivity is based on experimental evidence on cell level connectivity. Detailed reports of analysing this data as well as the computational modelling of inter-cellular connectivity are provided by publications: (Li et al., 2007;Borisyuk et al., 2011;Borisyuk et al., 2014;Roberts et al., 2014).

Locomotor central pattern generator (CPG).
Here we consider the neural populations that generate locomotor behaviour. Fig. 2 shows the diagram of connections between CPG populations. We consider seven interacting populations on each body side. The upper rectangle (light blue dotted borders) contains the neuronal populations relevant to swimming (Roberts et al., 2014): excitatory descending interneurons (dINs), inhibitory ascending interneurons (aINs) and inhibitory commissural interneurons (cINs). The lower rectangle (red dotted borders) contains the populations related to struggling activity (Li et al., 2007): excitatory descending repetitive interneurons (dINrs), inhibitory ascending repetitive interneurons (aINrs) and excitatory commissural interneurons (eCINs). In addition, we consider a population of motor neurons, which receive multiple inputs from other CPG populations and are considered the output of the CPG.
Note that experiments suggest that there is one aIN population, containing some neurons that are more reliably spiking during swimming and some that are more active during struggling. To simplify the model we split the aINs into two populations: dINs in the swimming circuit and dINrs in the struggling circuit.
The populations of dIN and dINr neurons can be considered as decision-making populations: they integrate input from the sensory pathways and select a corresponding action from the 6 repertoire. Swimming is characterised by the regular rhythmic activities of swimming populations with anti-phase oscillations between the equivalent populations on opposite body sides. The frequency of swimming is in the range 10-25 (Hz).

Model formulation: Sensory pathways
1. We assume that the touch skin (TS) population includes sensory RB, dla and dlc neurons. We combine these populations into one model population per body side. The average activity in these combined populations on the left and right sides is denoted by 2. Similarly, we assume that the TH population includes sensory tSt neurons and tINs, which are also combined together into one population per side. The average activity of this combined population on the left and right body sides is denoted by To describe activity dynamics of ten neural populations corresponding to the four sensory pathways and one integrating population (per side), we use a mean field approach based on the Wilson-Cowan equations. The system of ordinary differential equations (ODEs) describing the sensory pathway activities is: 10 10 4 4 2 3 10 10 10 10 10 5 9 3 2 1 1 9 9 9 9  Fig. 2 shows a schematic representation of the CPG populations and their connections. The swimming and struggling networks are shown by rectangles with dotted blue and red borders respectively. Each black arrow shows a directed connection from one population to another, where a sharp arrow end denotes an excitatory connection and a round end denotes an inhibitory one. Commissural connections travel through the vertical dashed line to the opposite body side. The violet arrows to the dIN and dINr populations indicate multiple inputs (either excitatory or inhibitory) from the sensory pathways on each body side that deliver signals to control action selection. The output from the CPG network is a pair of motor neuron populations on the left and right body side, which receive connections from all CPG populations on the same side. The activity in the motor neuron populations is considered to represent the tadpole's behaviour.

Model formulation: CPG
Thus, the CPG part of the model comprises three neuronal populations that are mostly active during swimming (dIN, aIN, cIN) on each side of the body. Their average activities we denoted by for the right side. The other three CPG populations (dINr, aINr, ecIN) on each side of the body are mostly active during struggling. We denote their average activities by The equations for the swimming related populations on the left and right body sides are: The equations for the struggling related populations on the left and right body sides are: Finally, the equations for the motor neuron populations are:  Note that we do not take into account the spatial distribution of neurons along the body but instead consider a set of localised populations of neurons. In the model formulation, we use a minimal set of populations with a minimal number of connections between them representing the strongest synaptic pathways.

Dynamics of the sensory pathway model
In this section, we analyse the dynamics of the populations that process sensory input. Processing of input for a particular modality occurs in a single population, and the equation governing this population's activity has the same form across all modalities. The same equation is also used to model the xIN populations, where input from the two touch pathways is integrated. The equation for a given population has the following form (extracted from equations (1)): Thus, the 1D non-linear equation (3) includes four parameters. Depending on the parameter values, there are three possible dynamical regimes: 1) a stable state with constant low activity; 2) a stable state with constant high activity; 3) bi-stability, where stable low and high activities coexist for the same parameter values. In the bi-stable regime, the system demonstrates hysteresis under parameter variation, and a short-term external input moves the system form one stable state to another. We use this bi-stable regime to model the dynamics of sensory pathways and decision-making populations of integrating xINs.
To find parameter values corresponding to the bi-stable regime we plot 2D bifurcation diagrams under variation of parameters: P a, keeping the other two parameters fixed. The bifurcation diagrams are shown in Fig. 3. merge. The area between these blue lines (shown by the horizontal line pattern in the rightmost bifurcation diagram) corresponds to a regime with three fixed points (two stable and one unstable). As the parameter values move towards the lines, the unstable (saddle) fixed point moves closer to one of the stable ones, and on the line the unstable and stable fixed points merge and disappear. Therefore, in the region outside of the lines there is a single stable fixed point, corresponding to either low or high constant activity depending on the parameter values.
To calculate the curve of saddle-node bifurcation on the co-dimension 2 bifurcation diagram, we used the software MATCONT (Dhooge et al., 2003) which is based on the idea of numerical continuation by parameters. In this particular case the saddle-node curve is described by two equations: equation 3 and another equation where the derivative of the right-hand side of equation (1) by the variable X equals to zero.
Generally the spiking activity of sensory neurons is noisy; we therefore include noise in the model of sensory populations (see formulas (1) and (3)): the external input is Here P is the constant part of input and  is a random Gaussian variable with zero mean and variance 2  . The value of the variance regulates the noise level; for example, the noise for pineal ganglion cells is substantially higher than the value of noise for other sensory pathways.
If we consider parameter values inside the region with three fixed points then if the variance of noise is small, e.g. 1 . 0   , the system stays at either the low or high level fixed point. If the variance is higher, e.g.
12 . 0   then the system will jump randomly between low and high levels of activity. To demonstrate this jumping we divide the interval of integration into a number of subintervals and on each subinterval with duration of 1 time unit, we randomly select a new value for the random Gaussian noise  . Fig. 4 the shows the resulting dynamics of population activity. This trajectory is typical for our model of sensory pathway activity.
If a sensory population receives external input (i.e. has a large value of ̅ ) then the population spends most of the time in the high (active) state. The population returns back to the mostly low (rest) state when the input decreases. Both touch skin (TS) and touch head (TH) populations can initiate swimming or struggling, and their activities are integrated by the xIN population. If at least one of HS or TH populations is active then the xIN population becomes active and sends excitation to the CPG populations to initiate locomotor behaviour.
The LD populations receive a sensory signal from the pineal eye. We assume that both populations receive a similar excitation and might become temporally active. The activity of LD populations on both sides converges onto the CPG populations and, if tadpole is not moving, this short-term signal can initiate swimming. If the tadpole is already swimming then the signal from LD populations leads to swimming acceleration, i.e. the period of regular oscillations in the swimming network becomes smaller. The detailed neurobiological mechanism of swimming acceleration is not yet known. To model this phenomenon we assume that the characteristic times ( 1 , … , 6 ) of the swimming populations are modulated by the LD activity level.
The PC population receives an input signal from the press head and cement gland sensors and becomes active for the duration of the sensory signal. The activity of the PC population inhibits dINs and stops swimming activity.
It is known from experimental study of tadpole behaviour that swimming can start spontaneously, apparently without any sensory stimulation (Jamieson & Roberts, 2000). To model this we modify the equations governing the xIN populations. We consider the external input to xIN populations as a function of time, which slowly increases when the animal is in the resting state. The slope of this increase is randomly and uniformly selected from a range, such that the resulting delays to spontaneous swimming start match biological evidence. When the activity in the xIN populations reaches a threshold then swimming starts, even without input from the TS or TH pathways. After swimming starts, the activity of the xIN populations returns to a low level and stays constant during swimming or struggling.
The neurobiological mechanism by which swimming can spontaneously stop is unknown (some hypotheses see at Dale & Giday, 1996;Zhang et al., 2014). To implement this behaviour we prescribe a maximum swimming duration. This duration is chosen randomly and uniformly from a suitable range. If the swimming is not disturbed by struggling then the swimming mode lasts until the end of the selected time interval.

Dynamics of CPG model
The CPG model can generate two main activity patterns: swimming and struggling. In a localised model the swimming pattern simply means the regular anti-phase oscillations of motor neurons (and other neuron types) on the left and right sides.
It is known from experiments that to initiate swimming a short touch of the skin on the trunk or head should be applied. This short touch provides a transient input (via the sensory pathway) to the dIN and dINr neurons on both sides of the body. Following the removal of sensory input, the CPG starts swimming and continues to work autonomously without additional signals from the sensory pathway. To model these dynamics we use a bi-stable regime, where a stable steady state (corresponding to low population activity / rest) coexists with a stable limit cycle (corresponding to anti-phase swimming oscillations) for the same parameter values. In this bi-stable regime, a short deviation of the system from the steady state leads to stable oscillatory activity.
To find appropriate parameter values for modelling swimming behaviour we used the results of bifurcation analysis of a Wilson-Cowan neural oscillator comprising two interacting populations (Borisyuk & Kirillov, 1992). A bi-stable region in parameter space has been found where a stable steady state and limit cycles (oscillations) co-exist, and we use parameter values from within this region. In a localised model, struggling is characterised by anti-phase bursting activity: the bursts arise in anti-phase on the left and right sides of the body and exist for the duration of the input signal. We model struggling behaviour as a regime of envelope oscillations, corresponding to a torus in the multidimensional phase space of the dynamical system. This dynamical regime is characterised by two frequencies of oscillations: a high frequency of intra-burst spiking activity and a low frequency burst "envelope". Motivation for using this type of dynamics comes from considering a population of more or less synchronously bursting neurons. The averaging of potentials of these bursting neurons will produce envelope oscillations with high and low frequencies. Example of such activity were considered in (Borisyuk 2002, Figure 2).
To find appropriate parameter values for modelling struggling behaviour we used the results of bifurcation analysis of two coupled Wilson-Cowan oscillators (Borisyuk et al., 1995). A regime of envelope oscillations was been found in this system and we use the corresponding parameter values for the model of struggling. Now we report results of CPG model simulations with the selected parameter values. We assume that without external influences the CPG model is in the rest state and the activity of all populations is near zero. To initiate swimming in the CPG model we apply a short external input to the dIN population on one side, e.g. on the left side. This approach is rather artificial; a more natural way of initiating swimming is via a relevant sensory pathway. However, it is common in neurobiology to activate a neural activity by direct stimulation and in our computer experiments we adopt this approach.
Experimental and modelling studies of swimming initiation follow trunk skin stimulation provide detailed neuronal mechanisms of this process (Roberts et al., 2014). The coarsegrained population model provides a simplified description of swimming initiation. We assume that stimulation of the dIN population on one side combined with a delayed stimulation of the dIN population on the opposite body side will initiate swimming. Fig. 5 shows simulation results where swimming is initiated by a short stimulus on one side. An external input 37 . 1 1  Q is applied to the dIN population on the left side from t=100 until t=250 ms (red horizontal bar in Fig. 5). To follow the experimental observations on appearance of activity on the opposite side, we apply an external input 3 . 1 4  Q to the dIN population on the right side with delay of 30 ms, i.e. from t=130 ms to t=250 ms (blue horizontal bar in Fig. 5). To break symmetry, we use slightly different stimulation durations and amplitudes on the two sides. After the stimulation period (i.e. after 250 ms) the external input to both dIN populations returns back to the initial value . This constant non-zero value of external input corresponds to the depolarization of dIN neurons due to their relatively high background NMDA activation (Roberts et al., 2014). Fig. 5 shows activity in the swimming populations (dINs, aINs, cINs) on the left (bottom traces) and right (top traces) in simulations with a short stimulus applied. Swimming activity starts on the side that is stimulated first, with the later stimulated side starting afterwards. The red and blue lines in Fig. 6 show activity in the motor neuron populations on the left and right sides respectively during the same simulation. It is clear that after a short transitional period with duration about 150ms the system generates regular anti-phase oscillations with period  T 50ms. Fig. 7 shows projections of two limit cycles to a 2D plane with like-wise populations on opposite sides. The larger limit cycle corresponds to antiphase oscillations in dIN populations on the left and right sides, while the smaller one shows anti-phase activity in cIN populations on the opposite sides. Fig 8 shows the activities of dIN and aIN populations on the same side. Fig. 8A shows a projection of the limit cycle to the 2D plane of dIN (horizontal) and aIN (vertical) axes. The shape of the limit cycle is typical for in-phase oscillations with a small temporal shift. Fig 8B shows oscillations of aINs (red) and dINs (blue). These two populations oscillate with a small phase shift of 10% of the period.
A simulation of struggling activity is shown in Fig. 9. Here, as before, we do not consider the sensory pathways but instead directly stimulate the dIN and dINr populations.  Fig. 9). In addition, we stimulate the dINr populations on both sides for 2000 1150   t (green stimulation bar). Because of this stimulation, modulated oscillations in the struggling network appear. The delay of 150ms in stimulation of the dINr populations in comparison with the dIN populations distinguishes swimming and struggling regimes, since struggling starts only in response to prolonged stimulation of both body sides. Without such delay every swimming initiation will also initiate activity in the dINr populations. Therefore, if stimulation of dINs is less than 150 ms long than dINr neurons are not stimulated and the swimming regime will be generated. In case the dIN stimulation is longer than 150 ms, dINr stimulation is applied and struggling behaviour is produced. Fig. 9 shows the activity in the motor neuron populations on the left (red) and right (blue) sides. Here the blue curve is shifted up by 0.07 units to avoid an overlap of graphs. The period of fast oscillations is about 5ms and the period of the slow (envelope) oscillations is about 130ms. These times are in line with experimental findings on struggling activity. Fig. 10 shows activities of dINr populations on the left (red) and right (blue) sides during the same struggling simulation. The blue graph is shifted up to avoid overlap. These graphs show envelope oscillations with slow anti-phase oscillations between the two body sides. Populations of aINrs and eCINs have similar dynamics during struggling (not shown). When stimulation is stopped (at time 2000 ms) the system returns to generating swimming activity. Fig. 11 demonstrates that struggling activity can also appear during swimming. To initiate struggling we first initiate swimming (as described above) and then we apply a long-term stimulation to dINs and dINrs for time 2000 1000   t , also as described above. Fig. 11 shows slow anti-phase envelope oscillations of the motor neurons on the left (red) and right (blue) sides. As before, struggling only persists during stimulation, with activity returning to swimming after the stimulation ends. The end of struggling therefore always leads to swimming, regardless of the state before the stimulation (rest or swimming).

Modelling of tadpole behaviour
In this section, we present simulation of the complete model with sensory pathways and CPG populations. Before running the simulation, we prepare a scenario that describes a sequence of external influences to the model's sensory inputs.
Here we describe one example scenario, which is summarized in Table 2. The scenario lasts from t=0 to t=7000ms. We assume that the tadpole is initially in the steady state (at rest); therefore, at t=0 all external inputs are initialized to zero. Event #1 is described as touching the skin on the left trunk during the time interval (200, 350) ms. Event #2 is light dimming in the time interval (1200, 1500) ms, etc. For each event in the scenario a relevant external input signal is generated. For example, for event #2, the input to TS sensory pathway changes for the duration of the event: the input to the TS population on the left side becomes 0.5 for 350 200   t and the input to the right-side TS population becomes 0.1 for 350 230   t . Event #2 causes swimming to start. Event #3 corresponds to light dimming, causing the input to the left and right LD populations to become non-zero for duration of event; as a result of this, swimming will be accelerated. However, it is possible that the swimming which was initiated by event #2 will stop spontaneously before the starting time of event #3. In this case the light dimming event will initiate swimming instead of accelerating it. Event #4 imitates a "predator attack", which is represented as a prolonged stimulation to both body sides. As a result of event #4, struggling behaviour occurs. Event #5 means a signal to pressure gland receptors, meaning the HP sensory pathway becomes active, causing inhibition of swimming behaviour and return to the steady state. Fig. 12 shows activity in the motor neuron populations vs time on the left (red) and right (blue) sides corresponding to the events in Table 1. Events are shown by vertical arrows. The duration of each event is shown by the horizontal bar under the time axis. These events represent a sequence of environmental signals to the tadpole's sensors, which modulate the dynamics of the corresponding sensory pathways. As a result of these influences from the sensory pathways, the CPG populations generate the relevant population dynamics. Motor neuron populations represent the output of the CPG and their activities demonstrate how the tadpole behaviour reflects the incoming environmental signals (events) which are specified in the scenario. The model dynamics are modulated by these signals to make decisions on appropriate locomotor actions and behaviour in response to environmental signals.

Discussion
In this paper, we have presented a computational model that is able to produce a set of behaviours mimicking the locomotor activities in a young tadpole (Roberts et al., 2010). These model behaviours arise in response to incoming signals from the "environment" according to a scenario consisting of a sequence of events. Thus, a complete meso-scale model of the tadpole nervous system has been developed. This model is able to process multiple sensory inputs and integrate them to make decisions about what locomotor response should be generated. The model is deeply rooted to data from many neurobiological experiments. For example, the mechanism of swimming initiation in response to skin touch has been studied and modelled at the level of spiking activity, and the neuronal populations participating in this process have been identified (Roberts et al., 2014;Buhl et al., 2015).
The model contains four sensory pathways as well as integrator and CPG populations. Although many experimental facts about the neurons in the sensory pathways and their electrophysiological properties are available, there are still many open questions and uncertainties, which require further investigation. Therefore, although we cannot derive a detailed model of spiking neurons for sensory pathways. We therefore use an approach that is based on a population model of Wilson-Cowan type. For consistency, the same approach has been used for modelling CPG activity.
Population level (mesoscopic-scale) models (also known as neural-mass models) are widely used in computational neuroscience. Several different approaches have been developed, but all consider neuronal activity in entire neuronal populations (Wilson & Cowan, 1972;Jancen et al., 1995;Marten et al., 2009;Liley et al., 2010;). In this approach, interacting populations influence each other by excitatory or inhibitory connections, which is usually represented by a sigmoid function (Marreiros et al., 2008). If the distribution of population's activity in space is considered then a generalization of neural-mass model known as a neural-field (mean-field) approach is applicable (Wilson & Cowan, 1973;Deco et al., 2008;Bressloff, 2012;Hlinka & Coombes, 2012). A natural extension to our model would be to consider spatial information, to attempt to reproduce the head-to-tail propagation of swimming activity that is seen in experiments, for example.
Decision making in our model is based on a slow increase of population activity toward a threshold. The activity of sensory populations is noisy and there is therefore variability in time of first threshold passage in the integrating population (Wang, 2002). This approach is in line with theoretical ideas on decision making in neuronal circuits (Kristan, 2008;Gold & Shalden 2007;Larsen & Bogatsz 2010;Marshal et al., 2012). Swimming in the model is based on a bi-stable regime. A short excitatory influence moves the system from a stable steady state to stable oscillations. Similarly, a short inhibition moves the system back to the resting steady state. Struggling modelling is characterised by bursting activities of participating neurons, therefore, it is appropriate to describe population activity during struggling by envelope oscillations with fast and slow frequencies.
Simulation of the model shows complex behaviour even in the case of a pre-defined scenario of input events. The reason is that the model includes spontaneous starting and stopping of swimming as well as random noise in the sensory populations. Therefore, repeated model simulations with the same parameter values according to the same scenario lead to different dynamics and behaviours. This makes the model behaviour less predictable and more similar to real tadpole movements in natural environments. Let us assume two robotic tadpoles each is controlled by the population model. Even if both robots are initiated at almost the same position and receive the same inputs, their behaviours might become significantly different after a short period of time. This fact means that the system is characterised by a complex and unpredictable behaviour, like a real tadpole. The population level model demonstrates a good correspondence between input events and behaviour. A further model development will include much more biologically realistic modelling of spiking neurons and reflect recent experimental findings regarding the functioning of sensory pathways and integrators (Buhl et al., 2015;Koutsikou et al., 2016). In addition, we would like to mention a possible implementation of the population model as a swimming robotic tadpole. An advantage of the population model over a detailed spiking neuron model is that the population model is much more suitable for robotic implementation since it requires much less computational power, making real-time simulation possible.

APPENDIX 1
Sigmoid functions of excitatory and inhibitory populations (indexed by e or i respectively) are: All other parameters and their value are given in Table 3. To break symmetries we randomise all parameters by adding white noise with variance equal to 0.01.
Parameter value choice is an important part of model development. We use three main approaches to this. (1) Parameter value is motivated by experimental measurements. This approach has been mostly used for parameter values of equations describing the activity dynamics of sensory pathways and references to experimental papers are given in Section 2.1. (2) Parameter value is selected according to theoretical (bifurcation theory) studies. We mostly use this approach to select parameter values for the swimming and struggling networks. Modelling of swimming is based on a bistable regime where a stable fixed point co-exists with a stable limit cycle. In (Borisyuk & Kirillov, 1993) we found this regime in a system describing a Wilson-Cowan oscillator (interacting excitatory and inhibitory populations) and starting from these basic parameter values we select values in the swimming network to provide the bistable regime. In a similar way, for modelling the struggling behaviour we use envelope oscillations with a high frequency related to spiking inside bursts (synchronous bursting on one body side) and a low frequency corresponding to the slow escaping movement. Of course, high frequency oscillations should be antiphase between the two body sides. In (Borisyuk et al., 2003) we found this regime of envelope oscillations in a system of two coupled Wilson-Cowan oscillators and we modify these basic parameter values to get anti-phase envelope oscillations in the struggling network. (3) Parameter value is selected as a result of multiple computational model simulations. This approach has been used in a few cases where neither inspiration from experiments nor from theory were available.       Table 1 summarises description of four sensory pathways. For each pathway the names of populations in that pathway are given, along with the names of corresponding dynamic variables in the model. These variables describe the dynamics of pathway populations.