https://orcid.org/0000-0002-1005-7671
Figures
Abstract
Paramecium is a large unicellular organism that swims in fresh water using cilia. When stimulated by various means (mechanically, chemically, optically, thermally), it often swims backward then turns and swims forward again in a new direction: this is called the avoiding reaction. This reaction is triggered by a calcium-based action potential. For this reason, several authors have called Paramecium the “swimming neuron”. Here we present an empirically constrained model of its action potential based on electrophysiology experiments on live immobilized paramecia, together with simultaneous measurement of ciliary beating using particle image velocimetry. Using these measurements and additional behavioral measurements of free swimming, we extend the electrophysiological model by coupling calcium concentration to kinematic parameters, turning it into a swimming model. In this way, we obtain a model of autonomously behaving Paramecium. Finally, we demonstrate how the modeled organism interacts with an environment, can follow gradients and display collective behavior. This work provides a modeling basis for investigating the physiological basis of autonomous behavior of Paramecium in ecological environments.
Author summary
Behavior depends on a complex interaction between a variety of physiological processes, the body and the environment. We propose to examine this complex interaction in an organism consisting of a single excitable and motile cell, Paramecium. The behavior of Paramecium is based on trial and error: when it encounters an undesirable situation, it backs up and changes direction. This avoiding reaction is triggered by an action potential. Here we developed an empirically constrained biophysical model of Paramecium’s action potential, which we then coupled to its kinematics. We then demonstrate the potential of this model in investigating various types of autonomous behavior, such as obstacle avoidance, gradient-following and collective behavior.
Citation: Elices I, Kulkarni A, Escoubet N, Pontani L-L, Prevost AM, Brette R (2023) An electrophysiological and kinematic model of Paramecium, the “swimming neuron”. PLoS Comput Biol 19(2): e1010899. https://doi.org/10.1371/journal.pcbi.1010899
Editor: Joseph Ayers, Northeastern University, UNITED STATES
Received: September 1, 2022; Accepted: January 26, 2023; Published: February 9, 2023
Copyright: © 2023 Elices et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Data Availability: Code for electrophysiology experiments, model fitting and figures can be found at https://github.com/romainbrette/Paramecium-model. Electrophysiology data and analyses, including particle image velocimetry can be found at https://doi.org/10.5281/zenodo.6074166. Behavioral data and code can be found at https://doi.org/10.5281/zenodo.6074480.
Funding: This work was supported by Agence Nationale de la Recherche (ANR-20-CE30-0025- 01 to RB, LLP, AMP and ANR-21-CE16-0013-02 to RB), Programme Investissements d’Avenir IHU FOReSIGHT (Grant ANR-18-IAHU-01 to RB), Fondation Pour l’Audition (Grant FPA RD-2017-2 to RB), CNRS (Défi Mécanobiologie, project PERCEE, to RB, LLP and AMP), and Sorbonne Université (Emergence, project NEUROSWIM to RB, LLP and AMP). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Competing interests: The authors have no competing interests.
Introduction
Behavior depends on a complex interaction between a variety of physiological processes, the body and the environment. This complexity makes it challenging to develop integrative models relating the different components. Thus, a strategy is to study model organisms that are experimentally accessible and structurally simpler than vertebrates. This strategy has been applied in particular to C. Elegans, with its 302 neurons and a known connectome [1–3]. However, electrophysiology is difficult owing to the small size of its neurons (about 3 μm), and developing empirically valid neuromechanical models of C. Elegans remains challenging. More recently, other model organisms have been introduced: Hydra, with a few thousand neurons and the advantage of being transparent [4,5], and jellyfish Aurelia aurita [6].
Here we propose to develop an integrative model of Paramecium tetraurelia [7], which is structurally much simpler than the abovementioned model organisms (Fig 1A), since it is a unicellular organism, while being large enough to perform intracellular electrophysiology (about 120 μm long [8]). Paramecium is a common ciliate, which swims in fresh water by beating its ~4000 cilia (Fig 1B) [8–10], and feeds on smaller microorganisms (bacteria, algae, yeast). It uses chemical signals to find food (Fig 1C and 1D), avoids obstacles thanks to mechanosensitivity (Fig 1E), displays collective behavior, adapts to changing environmental conditions and can even learn to respond to new stimuli [11].
A, Scanning electron microscopy image (left) [12], and bright field microscopy image (right) of Paramecium tetraurelia (scale bars: 10 μm). B, Typical helicoidal swim of Paramecium [13]. C, Accumulation of paramecia in a drop of acid [13]. D, Trajectory of a single cell in the acid drop, showing directional changes at the boundary [13]. E, Avoiding reaction against an obstacle, showing ciliary reversal followed by reorientation [14].
More than a century ago, Jennings described the basis of its behavior as “trial-and-error” [14]. Paramecium normally swims in a helicoidal fashion at about 500–1000 μm/s, but when it encounters something undesirable (obstacle, hot region, noxious substance), it produces an avoiding reaction (Fig 2A): it briefly swims backward, then turns and swims forward in a new direction. The avoiding reaction is triggered by a calcium-based graded action potential, which can be observed in an immobilized cell in response to a current pulse (Fig 2B). The calcium current is produced by L-type calcium channels located in the cilia [15], related to the CaV1 family found in neurons, heart and muscles of mammals [16]. Genes for many ionic channels have been found in the fully sequenced genome [17], and a number have been electrophysiologically identified [18]. Many signaling pathways of neurons have been found in Paramecium [19]. For this reason, Paramecium has been called a “swimming neuron” [20] and there is a vast amount of information about its electrophysiology, from studies done mainly in the 1960-80s [15,18]. However, there is no empirically based model of its action potential.
A, Typical spontaneous avoiding reaction: the ciliate swims backward, then turns and eventually resumes forward swimming, while spinning around its main axis during the entire movement. Images are separated by 150 ms, with intermediate shaded frames every 37 ms. The cell was placed in 20 mM NaCl and 0.3 mM CaCl2 to induce spontaneous avoiding reactions [21]. B, Intracellular recording of a voltage response (bottom right) to a square current pulse of amplitude 300 pA (top right) in an immobilized cell (left; A: anterior end; P: posterior end), showing a small action potential (in the standard extracellular solution, see Methods). The arrow points at a small upward inflexion due to the calcium current. Inset: Paramecium immobilized on a filter (background) with two electrodes. C, Velocity field of the fluid on a plane ~30 μm above the cell, calculated over the three shaded intervals shown in B. The blue arrow indicates mean velocity in the whole field, represented twice larger for clarity. The red arrows highlight the area neighboring the cell. C1, The fluid moves backward, which would make the cell swim forward. C2, The fluid moves forward. C3, The flow direction reverts on the posterior end, but not on the anterior right end, resulting in a swirling pattern.
We developed a biophysical model of Paramecium’s action potential, based on electrophysiological experiments on immobilized paramecia. We then augmented it by coupling ciliary calcium concentration with kinematic variables, with a phenomenological model constrained by simultaneous imaging of fluid motion induced by cilia beating and measurements of trajectories of freely swimming cells. In this way, we obtain a model of a swimming Paramecium, which exhibits spiraling movements and graded avoiding reactions. Finally, we show how the model can be used to investigate closed-loop behavior, including chemotaxis and collective behavior.
Results
A brief overview of Paramecium’s action potential
In this section, we recapitulate known facts about Paramecium’s electrophysiology, while illustrating our experimental techniques. To perform intracellular electrophysiology (see Methods, Electrophysiology), it is necessary to first immobilize the cell. To this end, we use a device we previously developed [22], which uses a transparent filter with holes smaller than the cell and a peristaltic pump. The pump draws the extracellular solution (4 mM KCl and 1 mM CaCl2) from an outlet below the filter, immobilizing the cells against the filter (Fig 2B, inset). Two high-resistance electrodes are then inserted into the cell and the pump is stopped. The cell is then held by the electrodes. One electrode is used to inject current, while the other is used to measure the membrane potential. Paramecium is a large cell (about 120 μm long and 35 μm wide for P. tetraurelia), which makes it isopotential [23,24].
Depolarization opens voltage-gated calcium channels located in the ~4000 cilia [8–10], similar to the L-type Cav1.2 family in mammals [16]. This can be noticed on Fig 2B as a small upward deflection before the peak of the membrane potential (arrow). This calcium current, denoted as ICa, activates rapidly (a few milliseconds), producing a current of up to about 4 nA [25]. Calcium entry then makes cilia reorient, which makes them beat forward.
To observe the beating direction, we use 1 μm tracer particles and estimate their velocity with particle image velocimetry (PIV, see Methods), on a plane about 30 μm above the cell (Fig 2C1–2C3, note that fluid velocity is much smaller than swimming velocity as it is measured above the cell [26]). Before the stimulation, the cilia beat backward to the right (Fig 2C1), which tends to make the cell move forward, with a spiraling movement over to the left [27,28], with the oral groove facing the spiral axis [13,29]. When calcium enters the cilia and calcium concentration reaches about 1 μM [30], cilia reorient and beat forward, slightly to the left (Fig 2C2), which makes the cell move backward. The calcium channels inactivate rapidly (a few milliseconds) through calcium-mediated inactivation: intraciliary calcium binds to calmodulin, which then inactivates the channels [18,31,32]. Calcium is then expelled by diffusion, buffering and pumps, in particular plasma membrane calcium pumps (PMCA) identified in the cilia [33,34]. After calcium has entered the cilia, voltage-gated K+ channels located in the basal membrane rapidly open, producing a delayed rectifier current IKd that damps the membrane potential (Fig 2B, just after the peak). Calcium also activates a smaller K+ current IK(Ca), which can be seen in the prolonged hyperpolarization after the stimulation [35,36].
After the stimulation, when calcium concentration has decreased below ~1 μM [30], cilia revert to backward beating (Fig 2C3). We notice a swirling pattern on Fig 2C3, which can be attributed to an asynchronous reversal of different groups of cilia. We will show how this relates to the change in swimming direction seen on Fig 2A.
We used current-clamp recordings and PIV measurements of fluid motion to build a model of the action potential together with electromotor coupling, summarized in Fig 3. We chose to use current-clamp rather than voltage-clamp recordings because good control is difficult to achieve with high resistance electrodes and several important processes are calcium-gated rather than voltage-gated, making voltage-clamp less relevant. Our modeling strategy is outlined in Fig 4. We start by determining passive properties (Fig 5) by using responses to hyperpolarized pulses, which do not trigger the calcium-based avoiding reaction. Then we isolate the delayed rectifier current IKd mechanically by deciliating the cells, thereby removing the voltage-gated calcium channels (Fig 6). We use these findings to build a complete model of the action potential with calcium dynamics coupled with ciliary reversal (as sketched in Fig 3), which we fit to electrophysiological responses of ciliated cells together with measurements of fluid motion (Fig 7). To couple this electrophysiological model with kinematics, we examine a simplified hydrodynamic model where cilia reverse with heterogeneous calcium sensitivity, as suggested by our PIV measurements and previous studies (Fig 8). This study allows us to infer the form of couplings between calcium concentration and kinematic variables, which we partially constrain with behavioral measurements of freely swimming cells (Fig 9). With these couplings, the electrophysiological model can then be used to simulate free swimming and behavioral responses to stimuli (Fig 9). Finally, we augment the model with elementary models of transduction and show simulations of behavior in structured environments, such as obstacle avoidance (Fig 10) and gradient following (Fig 11).
There is more Ca2+ outside than inside, and more K+ inside than outside. Calcium enters through ciliary voltage-gated channels as a current ICa. It then quickly inactivates these channels, forming a negative feedback loop. Calcium activates motor proteins, triggering ciliary reversal, as well as a ciliary K+ channel, producing an outward K+ current IK(Ca). The motor activation results in calcium concentration [Ca2+]i modulating kinematic parameters. Calcium is then expelled, in particular by plasma membrane Ca2+ ATPases (PMCA). Depolarization opens voltage-gated K+ channels in the basilar membrane, presumably near the cilium, creating a current IKd (delayed rectifier). A linear leak current is also present. Finally, an inward rectifier current IKir opens at very hyperpolarized voltages, which has little impact on the avoiding reaction.
Each number indicates the corresponding figure. The left column (5–7) is the electrophysiological modeling. The right column (8–9) is the modeling of couplings between electrophysiology and kinematic variables. Put together, we obtain a model of a freely swimming Paramecium (9), which is then augmented with elementary sensory transduction models to yield model simulations in structured environments (10–11).
A, Top: voltage responses of one cell to negative current pulses (I = 0 to -4 nA in 300 pA increments; dashed lines: start and end of pulses), in the standard extracellular solution (4 mM KCl and 1 mM CaCl2). The arrow points at an inflexion due to the inward rectifier current IKir. Bottom: model responses fitted to the data, showing the inferred reversal potential of K+ (EK) and the half-activation voltage VKir of the inward rectifier current. B, Current-voltage relationship over all cells (mean ± standard deviation, measured at pulse end) in 4 mM KCl (grey) and 0 mM KCl (blue). Removing K+ from the extracellular solution largely suppresses the inward current. C, Activation curve of the inward rectifier current in the fitted models. The current activates below EK (EL is leak reversal potential). The solid curve is the activation function with median parameters, the dashed curve is the activation function of the cell shown in A. D, Fitted parameters over n = 28 cells, grouped in passive parameters and inward rectifier parameters.
A, Current-voltage relationship in deciliated cells, showing a strong delayed rectifier current for depolarized voltages, and the inward rectifier current for hyperpolarized voltages. B, Top: voltage responses of one cell to positive current pulses (I = 0 to 4 nA in 300 pA increments). Bottom: responses of the two-gate Boltzmann model fitted to the data, showing the inferred half-activation voltage of the delayed rectifier current (dashed). C, Activation and time constant of the delayed rectifier current as a function of voltage in fitted models, with median parameters (solid) and for the cell shown in B (dashed). D, Statistics of fitted parameters (n = 16).
A, Voltage responses of a cell (top) to two sets of current pulses (bottom), from 0 to 5 nA (in 300 pA increments) and from 0 to 300 pA (in 25 pA increments). B, Ciliary response to the same currents, measured as the cosine of the mean angle of the velocity field, relative to the anteroposterior axis. C, Close up of an action potential triggered by a 1.5 nA current pulse, with the model fit (dashed). The arrow points at an upward deflection due to the calcium current. D, Ionic current calculated by subtracting the estimated leak current from the capacitive current. The inward current (I<0, shaded) corresponds to the calcium current. Integrating this current yields a calcium entry corresponding to a 10 μM increase in intraciliary calcium concentration. E, Responses of the fitted model. E1, Voltage responses. E2, Ciliary responses. E3, Voltage-gated calcium current ICa (top, negative traces), delayed rectifier K+ current IKd (bottom) and calcium-activated K+ current IK(Ca) (top, positive traces) in the fitted model. Currents are shown with the electrophysiological convention, i.e., I<0 is inward. E4, Intraciliary calcium concentration in the fitted model. The dashed lines show the ciliary reversal threshold and the half-inactivation concentration. F, Ionic currents inferred by the model for the action potential shown in C. G, Statistics of fitted parameters (n = 18).
A, Direction of fluid motion during forward swimming (blue) and backward swimming (red), relative to the anteroposterior axis. Averages are shown by arrows. B, Example of helicoidal motion of Paramecium, with the oral groove facing the axis. Highlighted frames are spaced by 750 ms. C, The translational velocity vector v is oriented along the anteroposterior axis. The rotation vector ω is in the dorsoventral plane (including the oral groove), making an angle θ with the anteroposterior axis. D, Rotating movement of the cell at the end of avoiding reactions of increasing strength [27]. E, Organization of ciliary basal bodies on the oral (ventral) side [10]. The oral apparatus (oa) is in the center (R: right; L: left; A: anterior; P: posterior; as: anterior suture; ps: posterior suture; cy: cytoproct). F, Calculation of kinematic parameters v, θ and ω in a spherical model of radius 60 μm, during successive phases of the avoiding reaction. First column: cilia beat to the rear and right, producing an axisymmetric force field pushing the organism forward while spinning around its axis. Local force amplitude is adjusted for a velocity of 500 μm/s. Second column: cilia revert and now beat to the front and right, pushing the organism backward. Third column: anterior left cilia revert back to the initial direction while anterior right cilia still beat towards the front, and posterior cilia partially revert, beating to the right. Translational velocity is now 0 and the rotation axis tilts to about 34°. Spinning speed ω also increases by a factor four. Fourth column: all cilia revert back to the initial beating direction. G, Measurement of fluid velocity in a sample cell beyond the anterior end (top) and beyond the posterior end (bottom), in response to 100 ms positive current pulses (1–5 nA), relative to the anteroposterior axis. H, Over n = 9 cells, the direction of posterior motion reverts back about 30 ms after anterior fluid motion (dashed line: linear regression; solid line: identity). Reversal duration is calculated as the time when cos(α) crosses 0, relative to the pulse end time.
A, Velocity v as a function of intraciliary calcium concentration [Ca2+] in the model. B, Angle θ of the rotation axis as a function of [Ca2+] in the model. C, Spinning speed ω as a function of [Ca2+] in the model. Angle and spinning speed increase at intermediate Ca2+ concentration, as implied by the spherical model in Fig 8F and 8D, Simulated model trajectory with three 2 ms current pulse stimulations of increasing amplitude. Images are shown at 400 ms intervals. Without stimulation, the organism swims in spiral, with the oral groove facing the spiral axis. A very small stimulation deviates the trajectory. Stronger stimulations produce avoiding reactions, with backward swimming and turning. E, Example of an observed Paramecium trajectory showing a directional change without backward swimming (right), followed by a full avoiding reaction (left). Images are shown at 400 ms intervals, starting on the right. F, Backward swimming duration (F1) and reorientation angle (F2) as a function of current amplitude for 2 ms pulses. Red and black curves show results for the same model but different initial positions of the oral groove, differing by a quarter of a cycle. G, Backward swimming duration (G1) and reorientation angle (G2) as a function of current pulse duration T with 100 pA amplitude. H, Reorientation angle vs. backward swimming duration in n = 1138 spontaneous avoiding reactions of Paramecium, showing a positive correlation (linear regression r = 0.2, p ≈ 10−11). About 15% of data points are not represented (larger angle or duration). Colors represent contour lines of the distribution.
A, Trajectory of the model doing several avoiding reactions against the stimulus. B, Membrane potential (top), stimulus current (middle) and intraciliary calcium concentration (bottom) at the first contact. Contact occurs at the boundary with the shaded region. Orange curves indicate backward swimming. Several weak avoiding reactions occur in succession. C, Trajectory of the model where sensory transduction has a 40 ms activation/deactivation time constant. In red, the stimulus is placed 300 μm further away. D, Same as B, for the black trajectory in C. The stimulus current lasts longer and peaks after the organism has started reacting, resulting in a stronger avoiding reaction.
A, Top: trajectories of 100 models swimming for 20 s in a torus with a depolarizing circular stimulus, modelled as in Fig 10C. The proportion of cells in the disk quickly decays (below). Membrane potential, stimulus current, and intraciliary calcium concentration are shown for the trajectory highlighted in red, which does an avoiding reaction against the disk after a number of spontaneous avoiding reactions. B, 100 model trajectories with a circular stimulus triggering an adapting hyperpolarizing current. Organisms tend to make avoiding reactions on the inner boundary of the disk. The proportion of cells in the disk increases over time. The highlighted trajectory enters the disk around t = 5 s with a large hyperpolarization, then displays several avoiding reactions against the boundary of the disk before exiting. C, Paramecia swimming in a linear stimulus gradient, modelled as in B. The position of 200 cells starting at position x = 5 mm is displayed every 5 s. D, Collective behavior in model paramecia induced by respiration and chemosensitivity. CO2 produced by cells is displayed in shades of grey (normalized to the spatial peak), and diffusion is simulated. CO2 concentration represents an attracting stimulus modelled as in B and C.
Passive properties
We start by estimating the passive properties (resistance, capacitance, reversal potentials) with model fitting techniques (see Methods, Electrophysiological modeling and Model optimization), and we compare with previous measurements in the literature. To this end, we use voltage responses to hyperpolarizing current pulses (duration 100 ms, amplitude 0 to -4 nA in 300 pA increments; Fig 5A, top). Such stimuli are known to trigger different voltage-gated currents: a small inactivating calcium current [37–39], and a strong inward rectifier K+ current IKir, i.e., which lets current pass mainly in the inward direction, below the reversal potential of K+ [40,41]. None of these currents are involved in the avoiding reaction and the action potential, and therefore they will not be included in the final model. In particular, the small calcium current is related to the escape reaction, an increase in swimming speed triggered by hyperpolarizing stimuli [42,43]. However, the inward rectifier current will allow us to estimate EK, the reversal potential of K+.
When the pulse intensity is strong, an inward current can be seen to activate after ~15 ms. The inward rectifier current has the property to activate mainly below EK [44]. This can be seen by removing K+ from the extracellular solution (making EK = −∞). With 4 mM extracellular KCl, hyperpolarization below about -60 mV activates a strong inward current, which is largely suppressed in 0 mM KCl (Fig 5B). After the pulse, the K+ current switches from inward to outward as it passes EK. This results in a change in decay speed, which is noticeable at about -60 mV in Fig 5A (top, arrow). We use this property to estimate EK.
To this end, we fit a biophysical model consisting of a linear leak current and an inward rectifier current IKir (Fig 5A, bottom) (see Methods, Electrophysiological modeling) using the Brian 2 simulator [45] with the model fitting toolbox [46], which applies differential evolution and gradient descent for least square estimation of model parameters.
We modeled the inward rectifier current as a non-inactivating current with Boltzmann activation, two gates and a fixed time constant: (Eqs (4) and (5); gKir is maximal conductance and nKir is the gating variable). Fig 5C shows the activation curve
of the cell shown in Fig 5A (dashed) and the curve with median parameters (solid), which confirms that the current activates essentially below EK.
We find that the resting potential V0 is -24.5 mV ± 10.6 mV (n = 28; median -22.5 mV; Fig 5D). Oertel et al. [25] previously reported about -23 mV with a slightly different extracellular solution. Capacitance C is 289 ± 75 pF. By comparison, P. caudatum, which is larger, has a capacitance of about 700 pF. Since P. caudatum is 200 μm long and 46 μm wide [47] and P. tetraurelia is 115 μm long and 34 μm wide [8], a simple scaling would predict a capacitance (200 x 46) / (115 x 34) ≈ 2.35 times smaller for P. tetraurelia, i.e. about 300 pF, which is consistent with our estimates.
Resistance at rest R (including the contribution from IKir) is 126 ± 62 MΩ. Finally, we find EK = -48 ± 10 mV, corresponding to an intracellular K+ concentration [K+]i = 29 ± 11 mM (based on the Nernst equation and given that [K+]o = 4 mM). This is consistent with estimates in the literature obtained with various methods, varying between 18 and 34 mM [40,48–50].
We briefly describe the parameter estimation results for IKir, even though these will not be further used, as only EK plays a role in depolarized responses. The estimates for total conductance gKir and half-activation voltage VKir are variable across cells, presumably because these parameters are not well constrained by the data (gKir and VKir cannot be estimated independently in the voltage region where channels are mostly closed). Nonetheless, the results confirm that IKir activates essentially below EK. Activation slope (kKir = 32 ± 9 mV) and time constant (τKir = 16 ± 4 ms) are better constrained. With the estimated parameters, the inward rectifier current contributes about 14% of the resting conductance (median; 16 ± 14%).
A model of the deciliated cell
Next, we analyzed the delayed rectifier current IKd responsible for repolarization. The ciliary calcium channel can be pharmacologically blocked with W-7, but this drug is toxic to Paramecium [51]. Instead, we isolated the delayed rectifier current mechanically by removing the cilia with ethanol [47,52] (see Methods, Deciliation). This procedure does not kill the cell, and cilia grow back after a few hours. It removes the voltage-gated calcium channels, which are located in the cilia [16,47], and thereby also removes the calcium-activated K+ current IK(Ca). In addition, it is no longer necessary to use the immobilization device (Fig 6A, inset). As can be seen on Fig 3, this procedure should leave only the delayed rectifier current IKd, activated upon depolarization, and the inward rectifier current IKir, activated upon hyperpolarization. These two currents are indeed observed in electrophysiological measurements (Fig 6A).
We fitted a Boltzmann model of the delayed rectifier current (see Methods, Electrophysiological modeling, Eqs (6), (7), (12) and (13)), IKd = gKdn2(EK − V) (gKd is the maximal conductance, n is the gating variable), to responses to 100 ms depolarizing current pulses (0 to 4 nA in 300 pA increments). This model turned out to fit the data as well as a Hodgkin-Huxley model, but with fewer parameters (see Methods, Model optimization). Fig 6B (bottom) shows responses of this model, and Fig 6C shows the activation curve and voltage-dependent time constant with median parameters and those for the cell shown in Fig 6B, with detailed statistics in Fig 6D. The delayed rectifier current activates at a median value of VKd ≈ 21 mV (30 ± 23 mV) with a slope kKd ≈ 7 mV (9 ± 6 mV). The time constant peaks at ~4.1 ms (4 ± 1.3 ms) at a voltage Vτ ≈ 23 mV (26 ± 48 mV), with a slope kτ ≈ 12 mV (14 ± 8 mV).
Based on these results, we further simplified the model by enforcing VKd = Vτ and kτ = 2kKd. This simplification slightly increases the fit error (1.82 vs. 1.8 mV; p = 0.009, two-tailed Wilcoxon test), but reduces the number of parameters. We used this simplified model in the full model of ciliated cells (leaving its parameters unconstrained).
The action potential
We now build a model of the action potential of ciliated cells, coupled with cilia reversal (Fig 7A). Cilia revert with very small depolarizations, of just a few mV [42]. For this reason, we used two sets of pulses, large pulses from 0 to 5 nA in 300 pA increments, and small pulses from -100 to 500 pA in 25 pA increments (Fig 7A, bottom; only positive currents are shown). Tip potentials could fluctuate between these two sets, therefore we aligned the traces to the median resting potential of -22 mV, and we fixed EK at its median value of -48 mV. Simultaneously, we seeded the extracellular medium with 1 μm tracer particles and imaged their motion at a frame rate of 30 Hz. Particle image velocimetry was then used to calculate the fluid velocity field, giving an indication of the direction of ciliary beating, as illustrated in Fig 2C (See Methods, Particle image velocimetry). Fig 7B shows the cosine of the mean angle α of the velocity field during stimulation, relative to the cell’s anteroposterior axis: 1 means that particles flow towards the posterior end, i.e., the cell is trying to swim forward; -1 means that the cell is trying to swim backward. Thus, cilia revert for a duration longer than the pulse (120–200 ms in this cell), graded with pulse intensity.
In the electrophysiological responses, we notice a small upward deflection after stimulation as shown in Fig 7C (arrow). This deflection is due to an inward current, the Ca2+ current. This current can be estimated by subtracting the estimated leak current from the capacitive current C.dV/dt (passive properties estimated by model fitting, see below). With a pulse of intensity I = 1.5 nA, we find that the inward part of that current peaks at about 2 nA (Fig 7D). Note that in this Figure, currents are shown with the electrophysiological convention so as to compare with previous voltage clamp studies [18], i.e., I<0 is inward; this is the opposite of the modeling convention. The peak amplitude shown in Fig 7D is an underestimation of the calcium current since part of the inward current may be masked by the K+ current, but it is comparable to previous estimations in voltage-clamp [25]. The calcium current is known to activate and inactivate quickly, within a few ms [25,31], as can be seen on Fig 7D.
We also observe small oscillations (Fig 7A), due to the interplay between Ca2+ and K+ currents, and a pronounced hyperpolarization after the pulse. This hyperpolarization is due to a calcium-activated K+ current IK(Ca). This current has been previously characterized electrophysiologically [35,36], as well as genetically and with immunochemistry [12,34].
Thus, we included the following currents in our model: a leak current IL, a voltage-gated calcium current ICa, with calcium-mediated inactivation, a delayed rectifier K+ current IKd, and a calcium-activated K+ current IK(Ca) (see Methods, Electrophysiological modeling). The calcium current ICa is produced by ciliary channels similar to L-type Cav1.2 channels [16]. We modeled it similarly to [32,53], but we allow for several inactivation binding sites (Eqs (14)–(16)). In addition, the current uses the Goldman-Hodgkin-Katz equation, which is more appropriate than the linear driving force (ECa—V) when intra- and extracellular concentrations are very different [54]. The calcium-activated K+ current IK(Ca) is simply modelled with a conductance increasing as a Hill function of calcium concentration (Eqs (17) and (18)).
In addition, the model must include calcium dynamics (Eq (19)). The decay of calcium concentration after the action potential may be due to a combination of processes, including diffusion towards the base, buffering (in particular by centrin), and pumps. We model this combination by a simple linear model of the calcium flux. However, this is not sufficient because there is a large calcium flux at rest through the voltage-gated calcium channels, which must also be expelled or buffered. Thus, we postulate that the resting calcium concentration, about 0.1 μM [55,56], is maintained by a pump operating near that concentration, which is consistent with properties of plasma membrane calcium pumps (PMCA), also present in the cilia [33,34]. We model this pump with Michaelis-Menten kinetics.
Finally, we couple calcium concentration with ciliary beating angle by a Hill function (Eq (20)).
We then fitted this complete model simultaneously to electrophysiological and motor responses to 100 ms current pulses (n = 18), while ensuring that the resting calcium concentration was 0.1 μM (Fig 7E). Thus, calcium concentration is not directly measured, but indirectly constrained by several processes: inactivation of ICa, activation of IK(Ca), and ciliary reversal.
Fig 7E1 and 7E2 show the fits for the cell shown in Fig 7A and 7B. Consistently with previous voltage-clamp measurements [25], the inferred calcium current is transient and peaks at about ~2.5 nA (-2.5 nA with the electrophysiology convention, Fig 7E3). A residual current remains, so that calcium concentration remains high in the model during stimulation (Fig 7E4). This is consistent with the fact that ciliary reversal can last for many seconds when the membrane is depolarized [42]. With large currents, inferred calcium concentration raises to about 22 μM, similar to previous estimations.
The voltage-gated potassium current is delayed relative to the calcium current, and the calcium-activated K+ current raises more slowly and is only dominant during repolarization (Fig 7E3), with a maximum of 0.65 nA. This is consistent with previous studies of that current [36]. Fig 7F shows the three different currents during the action potential shown in Fig 7C. As previously argued, the calcium-activated K+ current has a small contribution to the early current [25].
Over the n = 18 cells, we find in the fitted models that the calcium current has half-activation voltage VCa = -1 mV (0 ± 6 mV), activation slope k = 4.3 mV (4.3 ± 1 mV) and time constant about 0.9 ms (0.9 ± 0.4 ms) (Fig 7G). Estimated conductance is not well constrained and often very large. This is presumably because the peak current is mainly determined by the inactivation properties, and therefore the conductance parameter is not well constrained. Half-inactivation occurs at about KCa = 3.7 μM (log10(KCa in M) = -5.4 ± 0.2). This is close to patch-clamp measurements on cardiac L-type calcium channels [57]. The fitted models have about 4 binding sites (4.4 ± 1.7), larger than previous models [32,53] (which have a single site but were not constrained by Paramecium data). Calcium decays in the model with a median time constant of 130 ms and the pump operating near rest has a median maximum rate of 87 μM/s.
The delayed rectifier current has similar fitted parameters as in deciliated cells (median kKd = 4.9 mV, τKd = 4.5 ms), except half-activation tends to be lower (median VKd = 4 mV)–this might be because the responses are essentially below VKd, in the unsaturated part of the activation curve. The calcium-activated K+ current IK(Ca) has low affinity (KKCa = 63 μM, log10(KKCa in M) = -4.2 ± 0.7). This is consistent with the observation that in voltage-clamp, this current keeps on increasing for at least one second [36]. There are about nKCa = 3 binding sites is (3.5 ± 2.3).
Finally, cilia revert when the inferred calcium concentration is about 2.4 μM (log10(Kmotor in M) = -5.9 ± 0.8). This is close to measurements with triton-permeabilized cells, reporting about 1 μM [30]. We note that this and other concentration parameters depend on the estimation of intraciliary volume, which is approximate.
Overall, parameters of the fitted models are compatible with known properties of the currents and of ciliary reversal.
Swimming and turning
We now examine how Paramecium swims and turns, before coupling the electrophysiological model with swimming motion.
Before stimulation, the flow produced by the cilia is directed towards the posterior end, about 11° to the right (Fig 8A, blue). This should produce a forward left spiraling movement, as documented from observations of free swimming [27,28]. During a pulse that triggers an action potential, the flow is directed towards the anterior end, about 9° to the right (Fig 8A, red). This would make the cell swim backward, also spiraling to the left. An example of this spiraling motion is shown on Fig 8B: the organism swims forward while spinning around a tilted axis. The oral groove faces the spiral axis [13,29], which means that the rotational velocity vector ω is tilted from the main axis towards the oral side by an angle θ (in the median plane; Fig 8C). In freely swimming paramecia, we found that θ ≈ 13° (±6.4°) and the rotation speed ‖ω‖ is about 1 cycle/s (1.03 ± 0.2 cycle/s) (see Methods, Behavioral measurements and Behavioral analysis).
How does Paramecium turn? A directional change can occur if the angle θ changes. According to Jennings [27], θ increases during the avoiding reaction, in relation with stimulus strength (Fig 8D). In the final model, we will directly couple the calcium concentration [Ca2+] to the three kinematic variables v, θ and ω. In order to postulate plausible calcium-kinematic couplings, and to understand the relationship between ciliary beating patterns and kinematic parameters, in particular θ, we now examine an idealized hydrodynamic model, consisting of a sphere of radius 60 μm, for which we can use analytical formula relating forces and motion (see Methods, Hydrodynamic model). Fig 8E shows the ciliature on the ventral (oral) side, which appears to be organized around the oral groove. In the model, we postulate that cilia on the left (L), right (R), anterior (A) and posterior (P) sides may beat in different directions. The fluid produces local forces opposite to the direction of ciliary beating, and at low Reynolds number, the total force and torque map linearly to the translational velocity vector v and the rotational velocity vector ω in the cell coordinate system [58].
At rest (low calcium concentration), cilia beat towards the rear, slightly to the right (Fig 8F, first column), so that the fluid produces a force towards the front, slightly to the left. If the direction of ciliary beating is identical everywhere in spherical coordinates (that is, in terms of the cardinal directions North/South/East/West, Fig 12), then the force field over the sphere is symmetrical with respect to the main axis. This makes both the total force and the total torque align with the main (antero-posterior) axis, and therefore v and ω are also aligned with that axis, that is, θ = 0. The organism then moves forward, with a spinning movement around the axis. We adjust the force so that the velocity is 500 μm/s, which makes the sphere spin at about 1.8 Hz. Upon stimulation, when calcium concentration is high, cilia revert and beat forward (Fig 8F, second column), making the organism move backward.
Thus, the organism cannot turn unless there is some asymmetry in the ciliary beating pattern. Machemer [59] and Párducz [60] observed that during the turning phase, anterior and posterior cilia beat in different directions; Jennings [27] observed that left and right anterior cilia beat in different directions, where “left” and “right” are relative to the oral groove. In Fig 8F (third column), we examine what happens if cilia beat in a swirling pattern around the oral groove: the left anterior cilia beating backward, the right anterior cilia beating forward, and the posterior cilia beating to the right. This corresponds to what would happen near the calcium concentration threshold for global ciliary reversal, if cilia revert back first in the left anterior part, then in the posterior part, then in the right anterior part. The swirling pattern suggests that the cell is going to turn around an axis tilted from the main axis, in the plane separating the left and right sides, and this can be confirmed analytically. In this configuration, the sphere does not move along the main axis (v = 0 mm/s), because the net force along that axis is null, but it turns along an axis tilted by θ ≈ 34° from the main axis. The sphere also spins about 4 times faster (ω/(2π) ≈ 7.2 Hz)–this is of course qualitative since Paramecium is not a sphere. It returns to moving forward when all cilia revert back (Fig 8F, fourth column).
This pattern of ciliary reversal is suggested on Fig 2C. However, given that the particle flow was measured on a plane ~30 μm above the cell, that the cell could take different shapes and that the position of the oral groove was often difficult to estimate, it was generally not possible to determine the precise pattern of ciliary reversal empirically. Nevertheless, it is possible to demonstrate that cilia revert asynchronously. We measured particle flow separately in two regions of the field, beyond the anterior end, and beyond the posterior end (Fig 8G). This was only possible for 9 cells, where video quality was sufficient (see Methods, Particle image velocimetry). We measured the mean angle α of the flow field during stimulation with current pulses between 1 nA and 5 nA; weak pulses were not included because particle density tended to be lower due to sedimentation (weak pulses were recorded after strong pulses).
Fig 8G shows in one cell that posterior cilia (bottom) revert back after anterior cilia (top). Thus, before stimulation, when calcium concentration is low, anterior and posterior cilia beat in the same direction. Beating direction is also spatially homogeneous during stimulation, when calcium concentration is high. However, after stimulation, anterior and posterior cilia beat in different directions for a short period. This was a reproducible finding across the 9 cells: the posterior side tends to revert back slightly later than the anterior side (Fig 8H) (p = 7.10−17, one-tailed Wilcoxon test), with a mean delay of 26 ms (s.d. 29 ms). This confirms Párducz’ observations [60], which were based on electron microscopy of cells fixed during the avoiding reaction.
Modeling the avoiding reaction
This idealized hydrodynamic model shows that an asynchronous control of cilia across the membrane by calcium can make the organism turn, by transiently increasing the rotation angle θ and the spinning speed ω. The exact mechanism by which calcium might differentially activate the cilia is beyond the scope of this paper (see Discussion). Here we simply postulate phenomenological couplings between calcium concentration and kinematic parameters, as suggested by our analysis (note that the hydrodynamic model is not explicitly included in our final model, as it was only used to explore the relation between ciliary patterns and kinematics). As illustrated in Fig 8C, we will assume that the translational velocity vector v is aligned with the main axis, so that it is fully parameterized by the velocity v, and that the rotation vector ω lies in the plane of the oral groove, so that it is parameterized by its angle θ relative to the main axis and the spinning speed ω. We assume that all three kinematic parameters (v, θ, ω) are functions of intraciliary calcium concentration [Ca2+] Fig 7A–7C) (See Methods, Electromotor coupling, Eqs (21)–(23)).
We model velocity as a Hill function of [Ca2+], with threshold equal to the reversal threshold KCa and n = 2 sites, linearly scaled to match the maximum positive and negative velocities measured empirically Fig 7A; (Eq (21)). In freely swimming paramecia, we observed that the median velocity was 472 μm/s (521 ± 285 μm/s) for forward swimming and 370 μm/s (411 ± 200 μm/s) for backward swimming. Thus, in the model, we simply set both forward and backward maximum velocity at v = ±500 μm/s. Fig 9A shows the resulting function for the cell shown in Fig 7A–7F.
We model both θ and ω as bell functions of [Ca2+] Fig 9B and 9C; Eqs (22) and (23); see Methods, Electromotor coupling), peaking when [Ca2+] is near the global ciliary reversal threshold Kmotor (2.4 μM), as suggested by our analysis of the spherical model (Fig 8F). The minimum angle is taken from measurements of trajectories of freely swimming paramecia (θ ≈ 13°). For the maximum angle, we choose θ = 90° to account for the strongest avoiding reactions observed by Jennings (Fig 8D), corresponding to a rotation normal to the main axis. The minimum spinning speed is based on measurements (ω ≈ 1 cycle/s), and the maximum spinning speed is set to 4 times the minimum, as in the spherical model (Fig 8D).
In this way, we obtain a model in which all kinematic variables are coupled to the electrophysiological model. We note that, in contrast with the electrophysiological part, this model is only loosely constrained by measurements (see Discussion). We can then calculate organism motion from these variables, and thereby simulate behavior in an environment (see Methods, Kinematics). In the following, the model of one particular cell is chosen for illustration, the same cell as in Fig 7.
We first examine the trajectory of a model stimulated by 2 ms currents of varying amplitude (Fig 9D) (see Methods, Behavioral scenarios). Without stimulation, the organism swims in a helicoidal path. With a small stimulation amplitude (0.3 nA), the organism changes direction without swimming backward. At larger amplitude (0.5 nA), the organism swims backward for a very short time then turns and swims forward. When the amplitude is increased (5 nA), backward swimming is more noticeable. Directional changes without backward swimming do occur in freely swimming paramecia, as illustrated on Fig 9E: the organism first changes direction without swimming backward, then does an avoiding reaction.
In more detail (Fig 9F1), we observe that the cell swims backward when current intensity exceeds a threshold (here 372 pA), then backward swimming duration tends to increase with intensity. The reorientation angle following the stimulation changes continuously with stimulation strength, but in a complex way (Fig 9F2). In particular, small stimulations can trigger large turns without backward swimming. In addition, the directional change depends on the initial position of the oral groove: the red and black curves of Fig 9F2 correspond to the same cell but an oral groove position (spinning angle) differing by a quarter of a cycle. This occurs because the observation plane is fixed while the organism spins around its main axis, and because movement is constrained in the plane.
The characteristics of the avoiding reaction also depend on stimulus duration (Fig 9G). If the pulse amplitude is fixed (I = 0.1 nA) and its duration is increased, then the duration of backward swimming increases (Fig 9G1), and reorientation angle tends to increase for large durations but is not monotonous near threshold (Fig 9G2).
When we examine spontaneous avoiding reactions of freely swimming paramecia, we find that both backward swimming duration and reorientation angle vary broadly (156 ± 81 ms and 114 ± 66°, respectively) (Fig 9H), and there is a small although highly significant correlation (linear regression, r = 0.2, p = 10−11). Thus, backward swimming duration and reorientation angle are variable and not deterministically related.
A closed-loop behavioral model of Paramecium
We now use the model in closed loop, to describe how the interaction between organism and environment may give rise to behavior (see Methods, Behavioral scenarios), as opposed to motor responses to predetermined input stimuli (open loop). The following scenarios are relatively abstract and are meant mainly as illustrations of the possible uses of this integrative model to investigate the relation between physiology and behavior.
First, we consider an organism swimming towards a generic object, which triggers a depolarizing current when in contact with the membrane (for example a chemical substance, or hot water) (Fig 10). Thus, we simply consider that the stimulus current is proportional to the surface area in contact with the stimulus (see Methods, Sensory transduction). In contrast with previous situations, the stimulus is not pre-determined but depends on behavior. As Dewey pointed out (1896), “the motor response determines the stimulus, just as truly as sensory stimulus determines the movement.”. When the organism touches the object, a current is triggered, which depolarizes the membrane (Fig 10A and 10B). As the cell swims into the object, the current increases until an action potential is triggered. The cell then swims backward, moves out from the object and the current stops. Thus, the sensory current is necessarily small and short, because the organism withdraws as soon as the current reaches threshold. This results in a small avoiding reaction, and the organism bumps again repetitively against the object until it finally escapes (S1 Movie).
Larger movements can be obtained if sensory transduction has slower kinetics Fig 10C and 10D and S2 Movie). Here, the sensory current follows the stimulation with a time constant of 40 ms, modelled with first order kinetics (to simplify, the spatial spread of channels is not modeled; (Eq (24)). In this case, the sensory current keeps on increasing (slightly) after the organism has started swimming backward and it lasts longer. This results in a larger avoiding reaction.
Although the model is deterministic, the directional change can be described as pseudo-random. If the object is moved 300 μm further (Fig 10C), then the organism escapes with a larger angle. This is because the organism spins while swimming, so that its oral groove takes a different position when it touches the object.
We now consider that the object is a disk in a square environment (Fig 11A and S3 Movie). To avoid boundary effects, we consider that the environment has the topology of a torus (paramecia escaping to the left reappear to the right). To account for spontaneous avoiding reactions (occurring at a rate of about 0.2 Hz in our behavioral measurements), we added a noisy current to the membrane equation. The paramecia are modeled as in Fig 10C and 10D, with slow transduction. Fig 11A shows 100 trajectories starting from random positions and simulated for 20 s. The proportion of paramecia inside the disk decreases from 19 to 7%. Thus, the disk acts as a repelling stimulus. Occasionally, a paramecium gives an avoiding reaction against the boundary, makes a large turn and swims backward into the disk. In this case, it continues being stimulated and swims backward through the disk until it escapes. This peculiar behavior might be avoided if currents of opposite polarity were triggered when stimulating the rear, as is the case for thermal stimuli [61] and some chemical substances [62]. Traces shown in Fig 11A show the membrane potential V, the stimulus current I and the calcium concentration [Ca2+] for the trajectory shown in red, where we can distinguish a number of spontaneous avoiding reactions and one avoiding reaction against the disk.
Next, we examine how an object can act as an attractive stimulus Fig 11B and S4 Movie). This behavior can be obtained if the cell responds to a hyperpolarizing stimulus with adaptation, so that a depolarization is triggered when the stimulus stops. A simple model that exhibits this behavior is one where the stimulus triggers currents through two pathways with different kinetics: a hyperpolarizing current with fast kinetics, and a depolarizing current of equal magnitude with slow kinetics (Eqs (25)–(27)). In this way, a transient hyperpolarizing current is triggered when the stimulus switches on, and a transient depolarizing current is triggered when the stimulus switches off. This is shown on Fig 11B. The proportion of paramecia inside the disk increases with time (from 15% to 45% in 20 s). It can be seen that paramecia tend to aggregate inside the disk, mostly near the boundary. This is due to the small avoiding reactions but also to the curvature of the disk.
Using the same model, we then place the paramecia in a linear stimulus gradient (Fig 11C and S5 Movie); the environment is toric in the transversal dimension (up and down boundaries are glued). We place all paramecia at the same initial position, with random orientations. The population then rapidly ascends the gradient. This occurs because a depolarizing current is produced when the stimulus decreases, triggering an avoiding reaction. This behavior has been observed in Paramecium with thermal gradients [63,64], and shares similarities with bacterial chemotaxis [65].
Finally, we show an example of collective behavior (Fig 11D and S6 Movie). Paramecium produces CO2 by respiration, which acidifies its medium, and it is attracted by weak acids [66–68]. As a result, it can form aggregates, for example around a source of food or at the bottom of a depression slide [68]. We simulate the production of CO2 by paramecia and its diffusion in the medium (a square with torus topology), together with sensitivity to CO2 modeled in the same way as in Fig 9B and 9C (see Methods, Collective behavior). CO2 concentration is represented in Fig 11D as grey shades, after normalization. In this simulation, paramecia progressively form aggregates.
Discussion
Summary
We have built an integrative model of Paramecium that combines electrophysiology and motility. The model is informed by previous experimental literature and constrained by specific electrophysiology and trajectory measurements. It consists of an empirically constrained biophysical model of the action potential, with a coupling to kinematic parameters, which is more phenomenological. It can be simulated as a model of autonomous behavior in various environments.
The electrophysiological model was built by model fitting to current clamp data, with calcium-dependent properties indirectly constrained by ciliary reversal data. This method recovered properties of individual currents compatible with previous measurements obtained by different means. For example, the calcium threshold for ciliary reversal was estimated to be ~2 μM, the same order of magnitude as measured by varying extracellular calcium concentration in Paramecium with permeabilized membranes [30]. This is notable since calcium was not measured but only inferred from electrophysiology (indirectly through the estimation of the calcium current by the fitting procedure). The fitting procedure also determined that the calcium-dependent K+ current is small during the action potential, as previously determined with voltage-clamp experiments [25], but dominant after stimulation. The magnitude and time scale of calcium currents estimated by fitting were also compatible with voltage-clamp measurements [25]. Quantitative fitting allowed us to estimate additionally the calcium inactivation threshold (~3 μM) and the number of inactivation sites (~4).
By measuring ciliary induced flows during action potentials, we found that ciliary reversal is not synchronous across the cell, confirming previous observations obtained with electron microscopy [60]. We showed with a simple hydrodynamic model that asynchronous ciliary reversal allows the organism to turn, namely if the order of ciliary reversal follows a swirling pattern around the oral groove. From these findings, we built a phenomenological model of the coupling between calcium concentration and the main kinetic parameters, constrained with trajectory measurements (speed, angle).
The integrated model shows helicoidal swimming with graded avoiding reactions, where backward duration swimming and reorientation angle increase with stimulus strength or duration. As observed in spontaneous behavior, the model can also slightly reorient without swimming backward.
Behavior of the autonomous model is more complex than stimulus-response experiments, because the relation between sensory stimulus and motor response is circular. In particular, we noticed that the interaction with an object (e.g. a chemical substance) critically depends on the properties of sensory transduction. For example, efficient avoidance of the object requires persistent stimulation, e.g. with slow sensory activation/deactivation. Furthermore, sensory adaptation to a hyperpolarizing stimulus makes the stimulus attractive. Another possibility, which is non-stimulus-specific, is that excitability adapts to voltage changes [69], for which there is some electrophysiological evidence in Paramecium [24]. This can allow the model to follow a stimulus gradient. Finally, collective behavior can arise if organisms are sensitive to a substance that they produce. In summary, relatively complex behavior can be generated by the interaction of this simple “swimming neuron” with its environment. We note that we have not included any specific empirical model of transduction, which would be required in any empirical investigation of behavior using our model.
Limitations
This work has many limitations. First, ionic currents were measured simultaneously rather than in isolation, although we could isolate IKd by deciliation. This was partly for technical reasons (one cannot measure IK(Ca) while blocking the calcium current), and partly to ensure a global fit of the entire model to the action potential. Nonetheless, current overlap may cause difficulties for model fitting. For this reason, we strived to choose the simplest models that captured the phenomenology.
A second limitation is that calcium was not directly measured. Instead, it was indirectly constrained by several observed phenomena: calcium-dependent inactivation of ICa, calcium-dependent activation of IK(Ca), and ciliary reversal. Calcium imaging has been performed previously in Paramecium by pressure injection of a calcium indicator, in other contexts [56,70,71]. However, it is technically challenging to perform quantitative time-resolved measurements of ciliary calcium, because the cilia represent a small fraction of the total volume (2–3%) and beat at about 20 Hz.
Related to this limitation, our estimates of calcium-dependent parameters, for example the ciliary reversal threshold, depend on an estimate of the effective ciliary volume. Compared to our estimate, based on electron microscopy measurements, this effective volume may be reduced by crowding, or increased by fast buffering. Changing this parameter results in proportional changes in concentration parameters. However, the fact that the fitted ciliary reversal threshold is close to the threshold measured on permeabilized Paramecium suggests that our estimate (also used in [25]) was reasonable.
Another limitation is we did not measure ciliary beating directly, but rather its effect on the fluid. This was motivated by the fact that we were interested primarily in the movement induced by ciliary beating, as well as by technical reasons. The constraints of measuring fluid motion while the cell is impaled with electrodes made it difficult to consistently achieve good spatial resolution over the entire recordings. In future work, high speed imaging of ciliary beating could be used to determine the spatial pattern of ciliary reversal with higher precision, although contractions of the cell may complicate the analysis.
Because of these technical limitations, our model of electromotor coupling was highly simplified, restricted to a phenomenological relation between calcium concentration and three kinetic parameters. It could also be that this relation is not instantaneous, involving more indirect pathways.
Finally, we considered only generic rather than biophysical models of sensory transduction, and we did not consider mechanical or hydrodynamic interactions with objects [72–74]. More generally, the behavioral repertoire of Paramecium includes other aspects that we did not attempt to model, such as the escape reaction [42,43], contractions [75], trichocyst discharge [76,77] and gravitaxis [78].
How Paramecium turns
Since our model directly couples calcium concentration to kinetic parameters, it is not tied to any specific hypothesis about the ciliary beating pattern. However, turning is only possible if ciliary reversal is asynchronous, leading to a strongly asymmetrical ciliary beating pattern, otherwise the action potential would only trigger a back-and-forth movement in the direction of the main axis. We have shown that one possibility, compatible with the observed movement, is that cilia beat in a swirling pattern around the oral groove.
It is known that there is some structural and molecular heterogeneity between cilia, in particular between locomotor and oral cilia [9]. Whereas basal bodies are regularly placed on the dorsal side, they are spatially arranged on the ventral side with a characteristic pattern around the oral groove (see Fig 8E). The beating frequency during helicoidal swimming is also spatially heterogeneous [79], and tail cilia are also known to be immobile [80]. Ciliary heterogeneity is likely to be more complex than a distinction between oral and locomotor cilia, because when Paramecium is cut in two pieces below the oral groove, both pieces can turn in a similar way [81]. Such heterogeneity is likely a general feature of motile microorganisms, some of which can exhibit complex gaits [82].
Machemer [59] and Párducz [60] described an asynchrony between anterior and posterior ciliary reversal, with anterior cilia returning to their initial beating direction before posterior cilia, which we confirmed with our PIV measurements. However, this is not sufficient to produce turning: if both anterior and posterior beating patterns are axisymmetric, any combination of them would still produce movement along the main axis. Jennings mentions that there is also a reversal asynchrony between the anterior left and anterior right cilia, such that all anterior cilia transiently beat towards the oral groove during the avoiding reaction [27]. More detailed investigation is necessary to clarify this question.
Physiologically, asynchronous ciliary reversal can be due to differential calcium sensitivity, that is, the calcium threshold for reversal might vary across cilia. This is the implicit assumption of our model. It could also be that there are differences in calcium entry or removal across cilia (e.g. a gradient of calcium channel expression, or calcium pumps, or calcium buffering molecules). Some studies suggest that cyclic nucleotides may also differentially regulate the reversal threshold [83,84].
Previous models
We are aware of two previous attempts to model Paramecium’s action potential, neither of which was based on quantitative measurements. Hook and Hildebrand [85] used a calcium channel model with instantaneous transitions, an ohmic current-voltage relation (instead of GHK), an inactivation state accessible only from the closed state, and no voltage-dependent K+ channel (only a model of IK(Ca), which is not the major K+ current). Kunita et al. [86] used a Hodgkin-Huxley type model with voltage-dependent inactivation of calcium channels, which is not the main inactivation mode of this channel [31], even though the phenomenon exists on a slow timescale [87]. The model included two calcium channels (fast and slow), for which there is no electrophysiological support, and their relative activation was given as a function of time after stimulus start (i.e., it is not modeled). The calcium-dependent K+ current was not included. Neither electrophysiological model was fitted to data, and neither was coupled to a kinematic model.
Perspectives
The model could be improved by addressing the technical limitations listed above. In particular, it would be most enlightening to measure calcium concentration in the cilia at high temporal and spatial resolution, although it might require new technical developments. Further investigations should be carried out to understand in detail how Paramecium turns: to measure the spatial pattern of ciliary responses and to determine how this heterogeneity is achieved physiologically.
We have addressed only the avoiding reaction of Paramecium. The modeling effort could be completed by addressing other behavioral aspects, such as the escape reaction (increased ciliary beating speed upon hyperpolarization), which involves distinct hyperpolarization-activated channels [7]. Paramecium is sensitive to many sensory modalities, including temperature, various chemical substances, mechanical stimulation, light. Thus, the model should be completed by models of sensory transduction, as well as of mechanical interaction with objects. This would allow us to use the model to investigate the physiological basis of behavior of Paramecium in complex environments.
Finally, this work opens the perspective of addressing complex autonomous behavior in ecological environments, including adaptation, learning and problem solving [7], with a systemic modeling approach.
Materials and methods
Paramecium culture and preparation
Cultures of Paramecium tetraurelia were obtained from Éric Meyer, Institut de Biologie, Ecole Normale Supérieure, Paris, France. For electrophysiological experiments (at Institut de la Vision), paramecia were co-cultured with Klebsiella pneumoniae, where each week 1 mL of culture was reinjected into 5 mL of Wheat Grass Powder (WGP) buffer supplemented with 1 μL of beta-sitosterol. Cultures were kept at room temperature (about 20°C). Cells were harvested in the early stationary growth phase, between 3 and 5 days after feeding them. To wash and concentrate cells for experiments, a droplet of culture (approximately 600 μL) was placed in a narrow neck volumetric flask before adding extracellular solution used for electrophysiology (see below). Due to negative gravitaxis [88], paramecia tend to accumulate at the top of the solution. Thus, after approximately 10 min, a concentrated population of cells were retrieved from the top of the flask and placed in a microcentrifuge tube for at least 3h for adaptation [50,89]. The tube was shaken before collecting cells to perform an experiment.
The culture method differed slightly for the behavioral measurements with freely swimming paramecia, because these cultures were done in another lab (Laboratoire Jean Perrin). Instead, before an experiment, bacteria were first grown in 5 mL of WGP for 24 h at 27°C, then paramecia were grown by adding 1 mL of Paramecium culture and 1 μL of beta-sitosterol to the bacterized WGP, for 48 h at 27°C in the dark. About 0.4 mL of cell suspension were then pipetted from the top of the culture tube into 4 mL of extracellular solution (see Electrophysiology), at least 20 minutes before an experiment.
Swimming pools
Freely swimming paramecia were imaged at room temperature (~25°C) in square pools of side length 30 mm and depth 340 μm. These were obtained using micromilling and molding techniques. A Plexiglas mold, consisting of a square trench, is first milled with a square end mill of diameter 1 mm using a CNC micro-milling machine (Minitech, Machinary Corp., USA). Then a liquid mixture of Poly-DiMethyl Siloxane (PDMS, Sylgard 184, Dow Corning, USA) and its crosslinker (10:1 mass ratio) is poured onto the Plexiglas mold. It is immediately placed in a vacuum chamber for at least 1 h to remove any air bubbles. Crosslinking of the mixture is then obtained by placing the whole in an oven at 65°C for at least 4 h. Finally, the resulting transparent elastomer pool is gently peeled off the mold and put on a microscope glass slide. Prior to any experiment, the pool is exposed to an oxygen plasma for about 1 min to render the PDMS surface hydrophilic and prevent the trapping of air bubbles.
Behavioral measurements
For all behavioral experiments, about 500 μL of the cell suspension is pipetted into the pool with a concentration of 300–600 cells/mL. Trajectories are imaged at 50 Hz with a CMOS camera (Blackfly S BFS-U3-51S5M-C, Flir, USA, 2448x2048 pixels2, 10 bits), acquired with its dedicated acquisition software (Spinview, Flir, USA). A high magnification variable zoom lens (MVL12X12Z, Thorlabs) is used and yields a pixel size of 3.81 μm. The pool is uniformly illuminated with a dark field configuration, by placing ~10 cm beneath the pool a square LED panel (EFFI-SBL, Effilux, France), on top of which a fully opaque mask is positioned, partially covering the LED panel (typically ¼ of its surface). The LED panel produces a red light (wavelength λ = 625 nm) to minimize phototaxis [90,91]. Movies of the swimming paramecia are 200 s long (see S7 Movie).
To limit hard drive space, images are stored without their background with lossless compression (TIFF format). The background image is computed by taking for each of its pixels the minimum pixel intensity over the first 100 frames. It is then subtracted from each frame, and pixels with an intensity value below a threshold (automatically computed with the “triangle method”, see e.g. [92]) are set to 0.
Trajectories are extracted with the open source tracking software FastTrack [93], and manually inspected for corrections. Briefly, the software fits an ellipse to the cell’s shape, and disambiguates front and rear based on the asymmetry of the pixel histogram along the main axis. Trajectories shorter than 1 s and sequences where the cell is immobile are discarded. Errors in front/rear identification are automatically corrected as follows: when the cell turns by more than 20° over two successive frames, it is considered an error and the angle is flipped.
Trajectories with circling motions are also discarded. To this end, we calculate the proportion of the trajectory where the cell turns clockwise versus anti-clockwise (for trajectories longer than 4 s). These proportions should be balanced (0.5) for helicoidal trajectories. The trajectory is eliminated if these proportions differ by more than 0.1 from the expectation, with manual confirmation. In total, there were n = 554 selected trajectories.
Behavioral analysis
Analysis of helicoidal trajectories.
We manually selected 20 trajectories presenting clear helicoidal motion in the focal plane from 2 experiments, totaling 121 s. In each helicoidal trajectory, cell orientation γ(t) varies periodically with period T. We fitted γ(t) to a sinusoidal signal. We found T = 1.02 ± 0.27 s (mean ± s.d.), corresponding to a spinning speed cycle/s or about 2π/s in radians. The amplitude was θ = 13 ± 6.4°, the angle relative to the spiral axis.
Analysis of avoiding reactions.
An avoiding reaction is defined as a portion of trajectory during which the cell swims backward. This backward swimming is detected when the instantaneous motion vector m and the orientation vector o (posterior to anterior) point to opposite directions, i.e., m ⋅ o < 0. Avoiding reactions consisting of a single pair of frames were discarded. The mean frequency of spontaneous avoiding reactions was calculated as the number of avoiding reactions across all trajectories, divided by the total duration, yielding 0.18 Hz.
Only reorientation events involved in planar avoiding reactions were selected, based on measurements of the eccentricity of the ellipse that best fits the shape of paramecia. Whenever this eccentricity went below 0.8, the event was discarded. The total reorientation angle was obtained by summing all successive instantaneous reorientation angles during the entire avoiding reaction.
In Fig 7G, the 2D probability density of reorientation angle and backward duration was calculated with Gaussian kernel density estimation.
Deciliation
Deciliated cells were obtained by adding 96% ethanol to a tube containing the previously washed and adapted cells in the extracellular solution up to a final concentration of 5% (v/v) [52]. Then the tube was shaken for 2 min and left to rest for 1 min. Deciliated cells were collected from the lower half of the tube since they no longer accumulate at the top of the solution. Cilia start to grow back after approximately 30 min. Thus, as described in [52], in some experiments we blocked cilia regrowth by adding 10 mM of colchicine to the extracellular solution.
Electrophysiology
The extracellular solution used in all experiments contains 1 mM CaCl2, 4 mM KCl and 1 mM Tris-HCl buffer with pH 7.2, except for Fig 5B, where there was no KCl (blue curve). Microelectrodes of ∼50 MΩ resistance were pulled using a micropipette puller (P-1000, Sutter Instrument) from standard wall borosilicate capillary glass with filament (o.d. 1 mm, i.d. 0.5 mm, Harvard Apparatus). They were filled with a 1 M KCl solution using a MicroFil non-metallic syringe needle (MF 34G-5, World Precision Instruments); a few recordings were done with 3 M KCl (no particular change was noticed).
We used an upright microscope (LNScope, Luigs & Newmann) with two objectives, a 20× air objective (SLMPLN Plan Achromat, Olympus) to locate cells, and a 40× water immersion objective (LUMPLFLN, Olympus) with DIC contrast enhancement for electrophysiology and imaging.
Paramecia were immobilized using the device described in [22]. Briefly, paramecia are immobilized against a transparent filter (Whatman Cyclopore polycarbonate membranes; diameter 25 mm, pore diameter 12 μm) thanks to a peristaltic pump (Gilson Minipulse 3) that circulates the fluid from below the filter to above the device. Two microelectrodes are then lowered into the cell, and the pump is stopped. The cell is then held in place by the electrodes.
Electrophysiology recordings were performed using an amplifier with capacitance neutralization (Axoclamp 2B and Axoclamp 900 A, Molecular Device) and an analog–digital acquisition board operating at a sampling frequency of 40 kHz (USB-6343, National Instruments). Custom Python programs (https://github.com/romainbrette/clampy) were used to control the acquisition board.
Membrane potential was recorded with the reading electrode while 100 ms current pulses of various amplitudes were injected through the second electrode, with at least 1 s between successive trials.
Particle image velocimetry
To measure the flows induced by cilia beating, the bath was seeded with 1 μm silica or polystyrene particles (~0.2 mM) after paramecia were immobilized and the pump was stopped. Because of sedimentation, particle density was typically higher at the beginning of the experiment. Images were recorded at 30 Hz with a high-sensitivity CCD camera (Lumenera Infinity 3-6UR) over a 1392×1392 pixels region of interest surrounding the cell (8 bits depth, pixel width 0.178 μm). Frames were synchronized with electrophysiology recordings using a digital trigger.
Frames were preprocessed by removing the background (average image) and band-pass filtering (subtraction of two Gaussian filters with standard deviation 1 μm and 1.3 μm). Consecutive frames were then analyzed with particle image velocimetry (PIV) using the OpenPIV Python package (https://github.com/OpenPIV/openpiv-python.git), which calculates the velocity field using image cross-correlation. We used 50 μm windows with 2/3 overlap.
In each frame, we calculated the mean angle of the velocity vector over the entire field, using circular mean (argument of the mean complex unit vector; occasional missed frames were discarded). We then subtracted the angle of the anteroposterior axis. The position of anterior and posterior ends was measured manually. As the two ends can be visually ambiguous, they were automatically corrected (by swapping) when the flow measured before stimulation was directed towards the anterior end (indicating backward swimming).
In Fig 8A, for each cell we averaged the mean angle over all currents and over the 300 ms before stimulus (blue) or over the second half of the stimulus (red), for positive currents (<5 nA).
We also calculated the mean angle in the anterior and posterior regions as indicated in Fig 8G. Each region is a half-plane orthogonal to the main axis, starting at one end. For this analysis, we selected n = 9 cells with high quality video recordings and clear cell positioning, indicated by an absence of missed frames and a pre-stimulus flow deviating by less than 45° from the main cell axis. The average was restricted to responses to large pulses (1 to 5 nA), because those were recorded before the small pulses and therefore had a higher density of particles (due to sedimentation).
Electrophysiological modeling
In this section, we describe the biophysical models. The parameter values are obtained by fitting the models to the data (section Model fitting).
Electrode model.
All recordings were done with two electrodes, an injecting electrode and a reading electrode. Because of the capacitance and resistance of the injecting electrode, the current injected in the cell is a low-pass filtered version of the command current [94]. To estimate this current, we model the injecting electrode as a simple RC circuit and estimate its parameters Re and τe from responses to small pulses, assuming passive cell responses:
(1)
(2)
(3)
where V is the membrane potential, assumed identical to the potential of the reading electrode, V2 is the potential of the injecting electrode, Ie is the current injected in the cell, and ΔV accounts for a difference in tip potentials. The membrane equation (first equation) is a rough linear model of the cell, but only parameters Re and τe are used subsequently, to estimate Ie from I according to the last two equations (ΔV has no impact on Ie and therefore can then be discarded). In the 29 ciliated cells analyzed for passive properties, we found Re = 121 ± 8 MΩ and τe = 1 ± 0.9 ms.
Ionic currents.
Paramecium electrophysiology is reviewed in [18] and updated in [7,95]. Paramecium in an isopotential cell [23,24]. Thus, we consider a single membrane equation:
where C is membrane capacitance, IL = gL(EL − V) is the leak current, IKd is the delayed rectifier K+ current responsible for repolarization, ICa is the ciliary voltage-dependent Ca2+ current, IK(Ca) is the calcium-activated K+ current, IKir is the inward rectifier K+ current and I is a stimulating current. We did not include a few other electrophysiologically identified currents that are less relevant for this study, namely: Na+ [96,97] and Mg2+ [98,99] currents (since our extracellular solution does not contain these two ion species), and hyperpolarization-activated calcium currents responsible for the escape reaction [37–39], which we did not model.
The inward rectifier current IKir is a K+ current activated by hyperpolarization, most strongly below EK [40]. It is modeled as follows:
(4)
(5)
where p = 1 or 2 (p = 2 in the final version). We made this simple modeling choice because this current was only used as a way to infer the reversal potential EK. In particular, we did not include inactivation [41]. We also tested a version of the model where the linear driving force (EK − V) is replaced by the Goldman-Hodgkin-Katz expression [54], but it made no significant difference in fitting results.
For the delayed rectifier current IKd, we tested several models of the type:
(6)
(7)
We tested two classes of models. The Hodgkin-Huxley model is:
(8)
(9)
(10)
(11)
where
For numerical stability (near x = 0), we use this special function in the code rather than the explicit expression.
The Boltzmann model is:
(12)
(13)
The voltage-gated calcium current ICa is a calcium-inactivated current located in the cilia [18,31]. The corresponding channels have been genetically identified; they are similar to the Cav1 mammalian family (L-type), with a putative calmodulin binding site [16]. We model it similarly to [53,100]:
(14)
where m is the activation gating variable, h is the inactivation gating variable, and fGHK(V) is the normalized current-voltage relation of the open channel, given by the Goldman-Hodgkin-Katz equation [54]. In the Hodgkin-Huxley model, this relation is linear. However, with very different intracellular and extracellular calcium concentrations, a better model is the Goldman-Hodgkin-Katz equation. Resting intracellular concentration is about 50–200 nM [55,56], and rises to an estimated 20 μM during an action potential [25]. In contrast, extracellular concentration is 1mM in our experiments. Thus, we neglect intracellular concentration, which yields:
where F is the Faraday constant, R is the gas constant, and T = 293 K is temperature (20°C). Here, extracellular concentration has been lumped into gCa, which is now homogeneous to a current, while fGHK(V) is unitless and has no free parameter.
The activation gating variable is governed by:
(15)
while the inactivation gating variable is a Hill function of intraciliary calcium concentration [Ca2+]:
(16)
This is similar to the model of [53], except that the number of sites nCa is allowed to be greater than 1, because we found that this was necessary to fit our data.
A calcium-activated K+ current has been identified by comparison with Pawn mutants lacking functional voltage-activated calcium currents [36]. The current is largely reduced by EGTA. Genomic analysis indicates the presence of both BK and SK channels, with SK channels immunochemically identified in the cilia [12,34]. We model the current as follows, with activation as a Hill function of [Ca2+]:
(17)
(18)
Calcium dynamics.
Resting intracellular calcium concentration [Ca2+]0 has been estimated between 50 and 200 nM [55,56]. We chose [Ca2+]0 = 100 nM. Calcium enters the cilia when calcium channels open. The concentration increase is spatially uniform along the cilium [101]. It then decreases by three mechanisms: buffering, pumps, and diffusion. Buffering can occur with a variety of calcium-binding proteins, an important one being centrin, located in the infraciliary lattice, at the base of cilia [102,103]. Plasma membrane calcium pumps (PMCA) have been identified in the basal membrane with low affinity, around 10−7 M [104], and also in the cilia [33,34]. Suppressing the ciliary PCMAs by RNA interference prolongs backward swimming, which means that they are involved in the removal of calcium after an action potential. In principle, calcium can also diffuse to the basal cytosol. However, this has not been observed [105]. This might be because of cilia volume compared to the cell, or because calcium is buffered at the base of cilia. Both phenomena can be modeled by diffusion to the cilium base, with fixed resting concentration at the boundary.
We lump these diverse mechanisms into two simple processes: a linear process, with rate proportional to ([Ca2+]-[Ca2+]0), modelling diffusion and low-affinity buffers, and a high-affinity process operating near rest, with rate given by a Hill function of [Ca2+], modelling PMCAs or high affinity buffers. This results in the following equation:
(19)
where v is the volume of cilia and F is the Faraday constant. It can be seen that the role of the high-affinity process in this model is to counteract the calcium flow at rest, namely J = Irest/Fvcilia, while the low-affinity process independently tunes the rate of calcium removal after an action potential.
There are 3000–4000 cilia in P. tetraurelia [8–10], with the upper estimates likely including oral cilia. Each cilium is 10–12 μm long [23,106]. Each cilium is 270 nm wide but the fiber bundle is 200 nm wide [9]. This yields a total volume between 950 and 2750 μm3. We used the previous estimate vcilia = 1700 μm3 from [25], which is compatible with these bounds, but the uncertainty is large. In addition, the effective volume might be smaller because of crowding, or larger because of fast buffering. In practice, an error in the estimation of ciliary volume will translate into an equivalent change in all calcium binding constants (as well as λ and J). Binding constants and volume cannot be determined independently, because the (inverse) volume effectively acts as a unit for those constants.
For stability, the numerical implementation of calcium-dependent equations used equivalent versions written as a function of
where [Ca]0 = 0.1 μM is the resting concentration. For example, inactivation is rewritten as:
This equivalent change of variables avoids numerical issues when [Ca2+] approaches 0. The calcium dynamics equation rewrites as follows:
Electromotor coupling.
Cilia reorient when intraciliary calcium concentration reaches about 1 μM [30]. We model the ciliary angle as a Hill function of [Ca2+]:
(20)
where angles are relative to the anteroposterior direction (α = 0 when cilia beat to the rear), and Kmotor is the reversal threshold.
Velocity is modeled as an affine transformation of a Hill function with coefficient n = 2, changing sign at Kmotor:
(21)
where vmax = 500 μm/s is maximum velocity (both backward and forward), according to our measurements.
The angle θ of the rotation axis and the spinning speed ω are modeled as bell functions peaking at Kmotor:
(22)
(23)
We set θmin = 13° based on our measurements, and θmax = 90°, to allow for planar rotations as illustrated in Fig 8D. We set ωmin = 2π/s (1 cycle/s) based on our measurements and ωmin = 8π/s (2 cycles/s), based on the doubling found in the spherical model (Fig 8F). Thus, no extra free parameter is introduced.
Sensory transduction.
Instantaneous transduction. In Fig 10A and 10B where a well-delimited object (half-plane or disc) acts as a stimulus, we first calculate the proportion of the cell surface in contact with the object. To this end, cell shape is determined by the formula proposed by Zhang et al. [107]:
where x is the position along the major axis and y the position along the minor axis, a = 120 μm is cell length, b = 35 μm is cell width and β = 0.15 is an asymmetry factor. We then simply calculate the intersection of cell and object shapes (as pixel images). The stimulus current is then I = I0p, where I0 is maximum current and p is the fraction of the cell image within the object.
Delayed transduction. In Fig 10B and 10A, we simply consider that the transduction current activates and deactivates with a time constant τI:
(24)
This simple model corresponds to channels with finite opening and closing rates (namely, opening rate α = τIs and closing rate β = τI(1 − s), where s is proportional to the stimulus I0p) but does not take into account the spatial recruitment of channels.
Adaptation. In Fig 11B, 11C and 11D, we consider that there are two pathways with opposite polarity, a fast pathway and a slower pathway:
(25)
(26)
(27)
where τfast = 40 ms and τslow = 200 ms. Thus, for a constant stimulus s, the stationary current is 0. In Fig 11B (disc stimulus), the stimulus is s = I0p. In Fig 11C and 11D, where the environment is spatially continuous, the stimulus is simply the value at the center of the cell.
Model optimization
Model parameters are estimated with the model fitting toolbox of the Brian simulator [45,46] (https://github.com/brian-team/brian2modelfitting). Briefly, the software performs least square optimization using a combination of global optimization algorithms (we used differential evolution) and gradient descent, where the gradient is calculated symbolically from the model equations. Optimization with multiple objectives is done by adding the errors associated to the different objectives. Compiled code is automatically produced by code generation. Each fitting procedure took up to a few hours, and fitting scripts were run in parallel on different cells, using a small cluster of 3 PC with 8 cores each.
After model fitting, cells were discarded if passive properties were abnormal, indicating a bad recording (C > 500 pF or R<30 MΩ or R>500 MΩ), or if EK>EL (which is biophysically impossible), indicating a fitting problem.
Electrode fitting. First, for each cell we estimated electrode resistance Re and time constant τe from responses to small 100 ms pulses, both hyperpolarizing and depolarizing (|I|<0.5 nA), using Eqs (1)–(3). The error criterion was the sum of quadratic errors on both electrode potentials, measured from 100 ms before to 100 ms after the pulse. The estimated parameters were then used in subsequent fits.
Hyperpolarized fits. For Fig 5, we fitted the models described above with IL and IKir, (Eqs (4) and (5)), with least square minimization of the error on the reading electrode potential, taken from pulse start to 50 ms after the pulse. The stimuli were 100 ms pulses with amplitude between -4 and 0 nA in 300 pA increments.
We fitted the model with p = 1 and with p = 2 for IKir. Using two gates (n2) gave better fits than using one (n = 40; median error 2.3 vs. 2.6 mV; p = 8.10−5, one-tailed Wilcoxon test). We also fitted the n2 model with the ohmic driving force replaced by a Goldman-Hodgkin-Katz model, but it did not yield any significant improvement (p = 0.27, two-tailed Wilcoxon test). Selection criteria (see above) were passed by n = 28 cells. Statistics of fitted parameters are shown in Table 1.
Deciliated fits. For Fig 6, we fitted the models described above with IL and IKd, with least square minimization of the error on the reading electrode potential, taken from pulse start to 50 ms after the pulse. The stimuli were 100 ms pulses with amplitude between 0 and 4 nA in 300 pA increments. Parameters EK and C were taken from the previous fit to hyperpolarized responses. To make sure the short onset is well captured, the time interval is split in two windows: the first 30 ms and the rest of the response, and each window is equally weighted (meaning that a data point in the first window contributes more than a data point in the second window).
We tested the Hodgkin-Huxley (HH) type model (Eqs (6)–(11)) and the Boltzmann model (Eqs (6), (7), (12) and (13)). The two models performed similarly (n = 21; median 1.78 mV vs. 1.69 mV; p = 0.59, two-tailed Wilcoxon test). The median number of gates was 1.13 in the HH model (1.1–1.8, 25–75% interval) and 1.7 (1.3–2.7) in the Boltzmann model. In both models, the minimal time constant was very small (median 0.1 vs. 0.6 ms). Therefore, we chose a model with 2 gates (n2) and a nearly null minimal time constant (100 μs for numerical reasons). A Boltzmann n2 model gave similar results to a HH n2 model (p = 0.96, two-tailed Wilcoxon test) and had fewer parameters, while performing similarly to the unconstrained model (p = 0.79, two-tailed Wilcoxon test). Therefore, we chose the Boltzmann n2 model. Selection criteria (see above) were passed by n = 16 cells. Statistics of fitted parameters are shown in Table 2.
Fitting results motivated us to further simplify the model by enforcing VKd = Vτ and kτ = 2kKd. This corresponds to a simple biophysical model where opening and closing rates are of the form e±V/k. This simplification slightly increases the fit error (1.82 vs. 1.8 mV; p = 0.009, two-tailed Wilcoxon test), but has the advantage of reducing the parameter set to a single kinetic parameter, the maximum time constant (median 3.2 ms).
Ciliated fits. For Fig 7, we fitted the complete model to electrophysiological and ciliary responses to two sets of 100 ms pulses, a set of pulses between 0 and 5 nA in 300 pA increments, and a set of pulses between -100 and 500 pA in 25 pA increments. The complete model consisted of IL, the simplified n2 Boltzmann model of IKd with aKd = 0.1 ms, and
(Eqs (6), (7), (12) and (13)), ICa (Eqs (14)–(16)), IK(Ca) (Eqs (17) and (18)), calcium dynamics (Eq (19)) and electromotor coupling (Eq (20)).
To deal with possible shifts in tip potential between the two sets, we aligned all traces to a resting potential of -22 mV (the median resting potential), and EK was fixed at -48 mV (the median estimated EK).
The optimization error combined an error on the reading electrode potential, an error on the ciliary angle (mean angle of the PIV analysis), and an error on resting calcium concentration. For the voltage error, the response was divided in two equally weighted intervals: from pulse start to pulse end, and from pulse end to 500 ms after the end (to capture the post-stimulus hyperpolarization). The angle error was defined as the quadratic error on the corresponding unit vectors, which is equivalent to , and applied on the interval from 100 ms before the pulse to 500 ms after it. Finally, we ensured that the resting [Ca2+] was 0.1 μM by inserting an error on [Ca2+] on the interval from 100 ms before the pulse to the start of the pulse (ECa = ([Ca2+] − 0.1 μM)2). This effectively ensures J = Irest/Fvcilia.
Selection criteria (see above) were passed by n = 18 cells. Statistics of fitted parameters are shown in Table 3.
Statistics
Fitting results obtained with different models were compared using the Wilcoxon test. Statistics are given as mean ± standard deviation. In box plots, the box shows the first and third quartile with the median value inside, and the whiskers are the minimum and maximum values excluding outliers, which are shown as diamonds and defined as those at a distance exceeding 1.5 times the interquartile range from the box.
Hydrodynamic model
In Fig 8F, we calculated the motion vectors from patterns of ciliary beating on a sphere of 60 μm radius. The velocity vector U and the rotation vector Ω are given by:
where F is the external force and L is the external torque [58]. The matrix is called the mobility matrix. For a sphere of radius r, the mobility matrix is diagonal:
We consider that each patch of membrane is subjected to a force from the fluid, in the direction opposite to the ciliary beating direction, and we calculate the total force and torque for different ciliary beating patterns.
We use spherical coordinates (θ, ϕ) (see Fig 12), where θ = 0 corresponds to the North pole, considered as the anterior end, and ϕ = 0 is the meridian corresponding to the oral groove. Thus, a surface element is:
The local force is tangent to the sphere and oriented with an angle α, where α = 0 is the direction of the meridian, pointing South. Cartesian coordinates are chosen so that the x axis is dorso-ventral, the y axis is oriented left to right, and the z axis is posterior to anterior (i.e., South to North pole).
With these conventions, the local force expressed in Cartesian coordinates is:
where A is the amplitude of the force per unit area and α is its angle (to obtain this result, start from the North pole, rotate along the y axis by θ, then rotate along the z axis by α). Therefore, the total force is:
The local torque is τ = r × F where r is a radius. In Cartesian coordinates, we obtain:
The reported velocity is v = ‖U‖, the angle of the rotation axis θrotation is calculated in the xz plane, and the spinning speed is in cycle/s. Local force amplitude is identical in all ciliary beating patterns, and chosen so as to obtain a forward velocity of 500 μm/s.
Three different patterns are represented in Fig 8F. In the forward pattern, the local angle is uniform: α = -170°, corresponding to a beating direction of 10°, downward to the right. In the backward pattern, the local angle is α = -10°, corresponding to a beating direction of 170°, upward to the right (obtained by up/down symmetry). In both cases, the local force is axisymmetrical with respect to the anteroposterior axis. It follows that both Ftot and τtot are aligned with the z axis (anteroposterior). This can be seen in the formulas above by integrating with respect to ϕ, which yields 0 for the x and y coordinates.
In the turning pattern, the left anterior quarter (ϕ ∈ [−π, 0], ) follows the forward pattern (α = -170°), while the right anterior quarter (ϕ ∈ [0, π],
) follows the backward pattern (α = -10°) and the posterior half (θ ∈ [π/2, π]) has local forces pointing to the left (α = -90°), meaning cilia beating to the right. The posterior half generates a rotating pattern around the main axis (by axisymmetry), without translational movement. Each anterior quarter generates forces and torques
With αL = π − αR, we then find that the y and z components of the total torque vanish. Added to the torque generated by the posterior part, we obtain a torque vector in the xz plane, which is the plane of the oral groove, separating the cell in left and right parts.
Kinematics
We consider that the organism is an object moving by rigid motion. The organism is characterized by a position vector x, and an orientation matrix R defining the rotation of the reference frame (frame of the organism), so that a point y on the reference frame is mapped to Ry in the observer frame. The reference frame is chosen as in the spherical model above, so that z>0 points towards the anterior end while the x axis the dorsoventral axis.
Translational velocity is assumed to be in the posterior-anterior direction only: v = [0, 0, v], so that
Let ω be the rotation vector in the reference frame. We assume it is tilted from the main axis by an angle θ in the xz plane:
Thus, kinematics is determined by three variables v, ω and θ, which are functions of calcium concentration, as detailed in Electrophysiological modeling (subsection Electromotor coupling).
Over a time dt, the organism rotates by Ω(dt) = I + [ω]. dt, where
is the infinitesimal rotation matrix, such that [ω]y = ω × y. Therefore, the orientation matrix changes as R(t + dt) = RΩ(dt), giving:
For numerical reasons, we use quaternions instead of matrices [108], implemented with the Python packages quaternion and pyquaternion. Orientation is then represented by a unit quaternion q, and kinematic equations translate to:
where
is the pure quaternion with imaginary part v, and
Confinement to a plane. We constrain the organism to move in a plane. This is done simply by rotating the orientation vector at each time step so that it lies in the plane. Concretely, we calculate the orientation vector in the observer frame:
Then we rotate the orientation vector around the axis u that is orthogonal to both the z axis and p: u = (0, 0, 1) × p, by an angle θ = sin−1(pz/‖p‖). The final orientation quaternion is then q′ = Q(u, θ)q, where Q(u, θ) is the rotation of axis u and angle θ.
Behavioral scenarios
In all simulations, the electrophysiological model fitted to the cell shown in Fig 7 is integrated with Euler method and a time step of 0.1 ms, using Brian 2 [45]. Kinematics were integrated with a time step of 1 ms for simulations with stimuli with sharp boundaries (Figs 7, 8, 8B and 9A) and 2 ms for simulations with spatially continuous stimuli (Fig 9C and 9D).
Avoiding reaction. In Fig 9D, models are simulated in a plane with 2 ms pulse currents of amplitude 0.3, 0.5 and 5 nA, at 1 second intervals. In Fig 7F, the stimulus is a 2 ms current pulse of 0.01 to 10 nA amplitude. In Fig 7G, the stimulus is a 100 pA pulse of duration 0 to 100 ms. The reorientation angle is calculated as the change in angle before and after the stimulus, averaged over 1 second (the spinning period).
Interaction with a stimulus. In Fig 10 and all subsequent figures, trajectories are constrained to a plane. In Fig 10A and 10B, the stimulus is a half-plane. It produces an instantaneous depolarizing current, proportional to the fraction of the cell shape inside the stimulus (see Sensory transduction), with maximum 5 nA. In Fig 10C and 10D, the stimulus additionally goes through a low-pass filter with time constant 40 ms, representing the activation/deactivation rate of the receptors.
Repelling and attracting discs. In Fig 11A and 11B, the stimulus is a disc of radius 1 mm within a 4 mm torus. 100 trajectories are simulated for 20 s, with random initial positions. A noisy current is added to the membrane equation so as to produce spontaneous action potentials at the observed rate of 0.2 Hz. It is modeled as an Ornstein-Uhlenbeck process:
with τn = 20 ms and σn = 9 pA. Physiologically, this corresponds to the random opening of K+ channels [109]. This noise is included in all subsequent simulations.
In Fig 11A, the stimulus produces a depolarizing current with a 40 ms time constant, as in Fig 8C and 8D. In Fig 9B, the stimulus produces an adapting hyperpolarizing current (see Sensory transduction; τfast = 40 ms and τslow = 200 ms), with a maximum amplitude of 1 nA.
Gradient following. In Fig 11C, the environment has the topology of a cylindrical surface, i.e., circular in the small dimension (500 μm) and linear in the long dimension (25 mm). The stimulus is a linear gradient of 100 pA/mm, with transduction modeled with adaptation as for the attracting disc (Fig 9B).
Collective behavior. In Fig 11D, cells produce CO2 by respiration, which then diffuses and acidifies the fluid. This is modelled by the diffusion equation:
where S is the transduction current triggered by CO2, α is the production rate and D is the diffusion coefficient of CO2. In water at 25°C, D ≈ 0.002 mm2/s but we accelerate it by a factor 5 to speed up the simulation. The production rate is α = 100 pA/s in a square pixel of width 20 μm.
Supporting information
S1 Movie. Interaction of the Paramecium model with a depolarizing stimulus, corresponding to Fig 10A.
https://doi.org/10.1371/journal.pcbi.1010899.s001
(MP4)
S2 Movie. Interaction of the Paramecium model with a depolarizing stimulus and slow kinetics, corresponding to Fig 10C.
https://doi.org/10.1371/journal.pcbi.1010899.s002
(MP4)
S3 Movie. Model trajectories on a torus with a depolarizing circular stimulus, corresponding to Fig 11A.
https://doi.org/10.1371/journal.pcbi.1010899.s003
(MP4)
S4 Movie. Model trajectories on a torus with an attracting circular stimulus, corresponding to Fig 11B.
https://doi.org/10.1371/journal.pcbi.1010899.s004
(MP4)
S5 Movie. Model trajectories on a linear gradient (with circular topology in the vertical direction), corresponding to Fig 11C.
https://doi.org/10.1371/journal.pcbi.1010899.s005
(MP4)
S6 Movie. Model trajectories on a torus with CO2 production (breathing), diffusion, and chemosensitivity, corresponding to Fig 11D.
https://doi.org/10.1371/journal.pcbi.1010899.s006
(MP4)
Acknowledgments
We thank Marcel Stimberg for technical assistance and Eric Meyer for providing us with specimens of P. tetraurelia.
References
- 1. Cohen N, Denham JE. Whole animal modeling: piecing together nematode locomotion. Curr Opin Syst Biol. 2019;13: 150–160.
- 2. Cohen N, Sanders T. Nematode locomotion: dissecting the neuronal–environmental loop. Curr Opin Neurobiol. 2014;25: 99–106. pmid:24709607
- 3. Schafer WR. The Worm Connectome: Back to the Future. Trends Neurosci. 2018;41: 763–765. pmid:30366562
- 4. Dupre C, Yuste R. Non-overlapping Neural Networks in Hydra vulgaris. Curr Biol. 2017;27: 1085–1097. pmid:28366745
- 5. Wang H, Swore JJ, Sharma S, Szymanski JR, Yuste R, Daniel TL, et al. From neuron to muscle to movement: a complete biomechanical model of Hydra contractile behaviors. bioRxiv. 2020; 2020.12.14.422784.
- 6. Pallasdies F, Goedeke S, Braun W, Memmesheimer R-M. From single neurons to behavior in the jellyfish Aurelia aurita. Calabrese RL, Satterlie R, Kanso E, editors. eLife. 2019;8: e50084. pmid:31868586
- 7. Brette R. Integrative Neuroscience of Paramecium, a “Swimming Neuron.” eNeuro. 2021;8: ENEURO.0018-21.2021. pmid:33952615
- 8. Nagel U, Machemer H. Physical and physiological components of the graviresponses of wild-type and mutant Paramecium Tetraurelia. J Exp Biol. 2000;203: 1059–1070. pmid:10683165
- 9. Aubusson-Fleury A, Cohen J, Lemullois M. Ciliary heterogeneity within a single cell: the Paramecium model. Methods Cell Biol. 2015;127: 457–485. pmid:25837404
- 10. Iftode F, Cohen J, Ruiz F, Rueda AT, Chen-Shan L, Adoutte A, et al. Development of surface pattern during division in Paramecium. I. Mapping of duplication and reorganization of cortical cytoskeletal structures in the wild type. Development. 1989;105: 191–211.
- 11. Hennessey TM, Rucker WB, McDiarmid CG. Classical conditioning in paramecia. Anim Learn Behav. 1979;7: 417–423.
- 12. Valentine MS, Rajendran A, Yano J, Weeraratne SD, Beisson J, Cohen J, et al. ParameciumBBS genes are key to presence of channels in Cilia. Cilia. 2012;1: 16. pmid:23351336
- 13. Jennings HS. Studies on reactions to stimuli in unicellular organisms. ii.—the mechanism of the motor reactions of paramecium. Am J Physiol-Leg Content. 1899;2: 311–341.
- 14.
Jennings. Behavior of the lower organisms. New York, The Columbia university press, The Macmillan company, agents; [etc., etc.]; 1906. http://archive.org/details/behavioroflowero00jenn
- 15. Eckert R. Bioelectric Control of Ciliary Activity. Science. 1972;176: 473–481. pmid:5032346
- 16. Lodh S, Yano J, Valentine MS, Van Houten JL. Voltage-gated calcium channels of Paramecium cilia. J Exp Biol. 2016;219: 3028–3038. pmid:27707864
- 17. Martinac B, Saimi Y, Kung C. Ion Channels in Microbes. Physiol Rev. 2008;88: 1449–1490. pmid:18923187
- 18. Eckert R, Brehm P. Ionic mechanisms of excitation in Paramecium. Annu Rev Biophys Bioeng. 1979;8: 353–383. pmid:383005
- 19. Plattner H, Verkhratsky A. The remembrance of the things past: Conserved signalling pathways link protozoa to mammalian nervous system. Cell Calcium. 2018;73: 25–39. pmid:29880195
- 20.
Kung C, Saimi Y. Chapter 3 Ca2+ Channels of Paramecium: A Multidisciplinary Study. In: Bronner F, Kleinzeller A, editors. Current Topics in Membranes and Transport. Academic Press; 1985. pp. 45–66.
- 21. Kung C. Genic mutants with altered system of excitation in Paramecium aurelia. II. Mutagenesis, screening and genetic analysis of the mutants. Genetics. 1971;69: 29–45. pmid:5133084
- 22. Kulkarni A, Elices I, Escoubet N, Pontani L-L, Prevost AM, Brette R. A simple device to immobilize protists for electrophysiology and microinjection. J Exp Biol. 2020;223. pmid:32409484
- 23. Eckert R, Naitoh Y. Passive electrical properties of Paramecium and problems of ciliary coordination. J Gen Physiol. 1970;55: 467–483. pmid:5435781
- 24. Satow Y, Kung C. Voltage Sensitive Ca-Channels and the Transient Inward Current in Paramecium Tetraurelia. J Exp Biol. 1979;78: 149–161.
- 25. Oertel D, Schein SJ, Kung C. Separation of membrane currents using a Paramecium mutant. Nature. 1977;268: 120–124. pmid:593305
- 26. Drescher K, Goldstein RE, Michel N, Polin M, Tuval I. Direct Measurement of the Flow Field around Swimming Microorganisms. Phys Rev Lett. 2010;105: 168101. pmid:21231017
- 27. Jennings HS. The behavior of paramecium. Additional features and general relations. J Comp Neurol Psychol. 1904;14: 441–510.
- 28. Machemer H. Ciliary Activity and the Origin of Metachrony in Paramecium: Effects of Increased Viscosity. J Exp Biol. 1972;57: 239–259. pmid:5075893
- 29. Bullington WE. A further study of spiraling in the ciliate Paramecium, with a note on morphology and taxonomy. J Exp Zool. 1930;56: 423–449.
- 30. Naitoh Y, Kaneko H. Reactivated Triton-Extracted Models of Paramecium: Modification of Ciliary Movement by Calcium Ions. Science. 1972;176: 523–524. pmid:5032354
- 31. Brehm P, Eckert R, Tillotson D. Calcium-mediated inactivation of calcium current in Paramecium. J Physiol. 1980;306: 193–203. pmid:6257894
- 32. Eckert R, Chad JE. Inactivation of Ca channels. Prog Biophys Mol Biol. 1984;44: 215–267. pmid:6095365
- 33. Yano J, Valentine MS, Van Houten JL. Novel Insights into the Development and Function of Cilia Using the Advantages of the Paramecium Cell and Its Many Cilia. Cells. 2015;4: 297–314. pmid:26230712
- 34. Yano J, Rajendran A, Valentine MS, Saha M, Ballif BA, Van Houten JL. Proteomic analysis of the cilia membrane of Paramecium tetraurelia. J Proteomics. 2013;78: 113–122. pmid:23146917
- 35. Saimi Y, Hinrichsen RD, Forte M, Kung C. Mutant analysis shows that the Ca2+-induced K+ current shuts off one type of excitation in Paramecium. Proc Natl Acad Sci. 1983;80: 5112–5116. pmid:6410401
- 36. Satow Y, Kung C. Ca-Induced K+-Outward Current in Paramecium Tetraurelia. J Exp Biol. 1980;88: 293–304. pmid:7452140
- 37. Nakaoka Y, Iwatsuki K. Hyperpolarization-activated inward current associated with the frequency increase in ciliary beating of Paramecium. J Comp Physiol A. 1992;170: 723–727.
- 38. Preston RR, Saimi Y, Kung C. Calcium current activated upon hyperpolarization of Paramecium tetraurelia. J Gen Physiol. 1992;100: 233–251. pmid:1328468
- 39. Preston RR, Saimi Y, Kung C. Calcium-dependent inactivation of the calcium current activated upon hyperpolarization of Paramecium tetraurelia. J Gen Physiol. 1992;100: 253–268. pmid:1328469
- 40. Oertel D, Schein SJ, Kung C. A potassium conductance activated by hyperpolarization in paramecium. J Membr Biol. 1978;43: 169–185. pmid:712815
- 41. Preston RR, Saimi Y, Kung C. Evidence for two K+ currents activated upon hyperpolarization ofParamecium tetraurelia. J Membr Biol. 1990;115: 41–50. pmid:2110594
- 42. Machemer H. Frequency and directional responses of cilia to membrane potential changes inParamecium. J Comp Physiol. 1974;92: 293–316.
- 43. Roesle E. Die Reaktion einiger Infusorien auf einzelne Induktionsschläge. Z Für Allg Physiol. 1903; 139–168.
- 44. Lu Z. Mechanism of rectification in inward-rectifier K+ channels. Annu Rev Physiol. 2004;66: 103–129. pmid:14977398
- 45. Stimberg M, Brette R, Goodman DF. Brian 2, an intuitive and efficient neural simulator. Skinner FK, editor. eLife. 2019;8: e47314. pmid:31429824
- 46. Teska A, Stimberg M, Brette R. brian2modelfitting. Zenodo; 2020.
- 47. Machemer H, Ogura A. Ionic conductances of membranes in ciliated and deciliated Paramecium. J Physiol. 1979;296: 49–60. pmid:529122
- 48.
Hansma HG. Biochemical studies on the behavioral mutants of Paramecium aurelia: ion fluxes and ciliary membrane proteins. PhD Thesis, University of California, Santa Barbara. 1974.
- 49. Ogura A, Machemer H. Distribution of mechanoreceptor channels in theParamecium surface membrane. J Comp Physiol. 1980;135: 233–242.
- 50. Oka T, Nakaoka Y, Oosawa F. Changes in Membrane Potential during Adaptation to External Potassium Ions in Paramecium Caudatum. J Exp Biol. 1986;126: 111–117.
- 51. Hennessey TM, Kung C. An anticalmodulin drug, W-7, inhibits the voltage-dependent calcium current in Paramecium caudatum. J Exp Biol. 1984;110: 169–181. pmid:6747535
- 52. Ogura A. Deciliation and Reciliation in Paramecium after Treatment with Ethanol. Cell Struct Funct. 1981;6: 43–50.
- 53. Standen NB, Stanfield PR. A Binding-Site Model for Calcium Channel Inactivation that Depends on Calcium Entry. Proc R Soc Lond B Biol Sci. 1982;217: 101–110. pmid:6131420
- 54. Hille B. Ion Channels of Excitable Membranes. Sinauer Associates; 2001.
- 55. Iwadate Y. Photolysis of caged calcium in cilia induces ciliary reversal in Paramecium caudatum. J Exp Biol. 2003;206: 1163–1170. pmid:12604576
- 56. Klauke N, Plattner H. Imaging of Ca2+ transients induced in Paramecium cells by a polyamine secretagogue. J Cell Sci. 1997;110: 975–983. pmid:9152023
- 57. Höfer GF, Hohenthanner K, Baumgartner W, Groschner K, Klugbauer N, Hofmann F, et al. Intracellular Ca2+ inactivates L-type Ca2+ channels with a Hill coefficient of approximately 1 and an inhibition constant of approximately 4 microM by reducing channel’s open probability. Biophys J. 1997;73: 1857–1865. pmid:9336181
- 58. Lauga E, Powers TR. The hydrodynamics of swimming microorganisms. Rep Prog Phys. 2009;72: 096601.
- 59. Machemer H. Regulation der Cilienmetachronie bei der “Fluchtreaktion” von Paramecium*. J Protozool. 1969;16: 764–771.
- 60.
Párducz B. Ciliary Movement and Coordination in Ciliates. In: Bourne GH, Danielli JF, editors. International Review of Cytology. Academic Press; 1967. pp. 91–128.
- 61. Tominaga T, Naitoh Y. Membrane Potential Responses to Thermal Stimulation and the Control of Thermoaccumulation in Paramecium Caudatum. J Exp Biol. 1992;164: 39–53. pmid:32539248
- 62. Oami K. Distribution of chemoreceptors to quinine on the cell surface of Paramecium caudatum. J Comp Physiol A. 1996;179: 345–352.
- 63. Mendelssohn M. Ueber den Thermotropismus einzelliger Organismen. Arch Für Gesamte Physiol Menschen Tiere. 1895;60: 1–27.
- 64.
Jennings HS. Behavior of the lower organisms. New York, The Columbia university press, The Macmillan company, agents; [etc., etc.]; 1906. http://archive.org/details/behavioroflowero00jenn
- 65. Berg HC. Bacterial behaviour. Nature. 1975;254: 389–392. pmid:1090851
- 66. Dryl S. Chemotaxis in Ciliate Protozoa. Behaviour of Micro-organisms. Springer US; 1973. pp. 16–30.
- 67. Houten JV. Two mechanisms of chemotaxis inParamecium. J Comp Physiol. 1978;127: 167–174.
- 68. Jennings HS. Studies on reactions to stimuli in unicellular organisms. iv.—laws of chemotaxis in paramecium. Am J Physiol-Leg Content. 1899;2: 355–379.
- 69. Platkiewicz J, Brette R. A threshold equation for action potential initiation. PLoS Comput Biol. 2010;6: e1000850. pmid:20628619
- 70. Iwadate Y, Katoh K, Asai H, Kikuyama M. Simultaneous recording of cytosolic Ca2+ levels inDidinium andParamecium during aDidinium attack onParamecium. Protoplasma. 1997;200: 117–127.
- 71. Iwadate Y, Kikuyama M. Contribution of Calcium Influx on Trichocyst Discharge in Paramecium caudatum. Zoolog Sci. 2001;18: 497–504.
- 72. Berke AP, Turner L, Berg HC, Lauga E. Hydrodynamic attraction of swimming microorganisms by surfaces. Phys Rev Lett. 2008;101: 038102. pmid:18764299
- 73. Jana S, Eddins A, Spoon C, Jung S. Somersault of Paramecium in extremely confined environments. Sci Rep. 2015;5. pmid:26286234
- 74. Ohmura T, Nishigami Y, Taniguchi A, Nonaka S, Manabe J, Ishikawa T, et al. Simple mechanosense and response of cilia motion reveal the intrinsic habits of ciliates. Proc Natl Acad Sci. 2018;115: 3231–3236. pmid:29531024
- 75. Nakaoka Y, Machemer H. Effects of cyclic nucleotides and intracellular Ca on voltage-activated ciliary beating in Paramecium. J Comp Physiol A. 1990;166: 401–406.
- 76. Hamel A, Fisch C, Combettes L, Dupuis-Williams P, Baroud CN. Transitions between three swimming gaits in Paramecium escape. Proc Natl Acad Sci U S A. 2011;108: 7290–7295. pmid:21464291
- 77. Knoll G, Haacke-Bell B, Plattner H. Local trichocyst exocytosis provides an efficient escape mechanism for Paramecium cells. Eur J Protistol. 1991;27: 381–385. pmid:23194850
- 78. Jensen P. Ueber den Geotropismus niederer Organismen. Arch Für Gesamte Physiol Menschen Tiere. 1893;53: 428–480.
- 79. Jung I, Powers TR, Valles JM. Evidence for two extremes of ciliary motor response in a single swimming microorganism. Biophys J. 2014;106: 106–113. pmid:24411242
- 80. Machemer H, Machemer-Röhnisch S. Mechanical and electric correlates of mechanoreceptor activation of the ciliated tail inParamecium. J Comp Physiol A. 1984;154: 273–278.
- 81. Jennings HS, Jamieson C. Studies on Reactions to Stimuli in Unicellular Organisms. X. The Movements and Reactions of Pieces of Ciliate Infusoria. Biol Bull. 1902;3: 225–234.
- 82. Wan KY. Synchrony and symmetry-breaking in active flagellar coordination. Philos Trans R Soc Lond B Biol Sci. 2020;375: 20190393. pmid:31884920
- 83. Noguchi M, Ogawa T, Taneyama T. Control of ciliary orientation through cAMP-dependent phosphorylation of axonemal proteins in Paramecium caudatum. Cell Motil. 2000;45: 263–271. pmid:10744859
- 84. Noguchi M, Nakamura Y, Okamoto K-I. Control of ciliary orientation in ciliated sheets from Paramecium–differential distribution of sensitivity to cyclic nucleotides. Cell Motil. 1991;20: 38–46.
- 85. Hook C, Hildebrand E. Excitation of paramecium. J Math Biol. 1979;8: 197–214.
- 86. Kunita I, Kuroda S, Ohki K, Nakagaki T. Attempts to retreat from a dead-ended long capillary by backward swimming in Paramecium. Front Microbiol. 2014;5. pmid:24966852
- 87. Hennessey TM, Kung C. Slow inactivation of the calcium current of Paramecium is dependent on voltage and not internal calcium. J Physiol. 1985;365: 165–179. pmid:2411920
- 88. Naitoh Y, Eckert R. Electrophysiology of Ciliate Protozoa. Exp in Physiol and Biochem. 1972. pp. 17–31.
- 89. Machemer-Röhnisch S, Machemer H. A Ca paradox: Electric and behavioural responses of Paramecium following changes in external ion concentration. Eur J Protistol. 1989;25: 45–59. pmid:23195786
- 90. Iwatsuki K, Naitoh Y. Photobehavior in a Colorless Paramecium. BioScience. 1983;33: 714–715.
- 91. Iwatsuki K, Naitoh Y. Photoresponses in colorlessParamecium. Experientia. 1982;38: 1453–1454.
- 92. Zack GW, Rogers WE, Latt SA. Automatic measurement of sister chromatid exchange frequency. J Histochem Cytochem. 1977;25: 741–753. pmid:70454
- 93. Gallois B, Candelier R. FastTrack: An open-source software for tracking varying numbers of deformable objects. PLOS Comput Biol. 2021;17: e1008697. pmid:33571205
- 94. Brette R, Piwkowska Z, Monier C, Rudolph-Lilith M, Fournier J, Levy M, et al. High-resolution intracellular recordings using a real-time computational model of the electrode. Neuron. 2008;59: 379–91. S0896-6273(08)00539-4 pmid:18701064
- 95. Valentine MS, Van Houten J. Ion Channels of Cilia: Paramecium as a Model. J Eukaryot Microbiol. 2022; e12884. pmid:34995386
- 96. Saimi Y. Calcium-dependent sodium currents inParamecium: Mutational manipulations and effects of hyper- and depolarization. J Membr Biol. 1986;92: 227–236.
- 97. Saimi Y, Ling KY. Calmodulin activation of calcium-dependent sodium channels in excised membrane patches of Paramecium. Science. 1990;249: 1441–1444. pmid:2169650
- 98. Preston RR. Transmembrane Mg2+ Currents and Intracellular Free Mg2+ Concentration in Paramecium tetraurelia. J Membr Biol. 1998;164: 11–24. pmid:9636240
- 99. Preston RR. A magnesium current in Paramecium. Science. 1990;250: 285–288. pmid:2218533
- 100. Chad J, Eckert R, Ewald D. Kinetics of calcium-dependent inactivation of calcium current in voltage-clamped neurones of Aplysia californica. J Physiol. 1984;347: 279–300. pmid:6323696
- 101. Iwadate Y, Suzaki T. Ciliary reorientation is evoked by a rise in calcium level over the entire cilium. Cell Motil. 2004;57: 197–206. pmid:14752804
- 102. Plattner H. Molecular aspects of calcium signalling at the crossroads of unikont and bikont eukaryote evolution—the ciliated protozoan Paramecium in focus. Cell Calcium. 2015;57: 174–185. pmid:25601027
- 103. Plattner H, Klauke N. Calcium in ciliated protozoa: sources, regulation, and calcium-regulated cell functions. Int Rev Cytol. 2001;201: 115–208. pmid:11057832
- 104. Wright MV, van Houten JL. Characterization of a putative Ca2+-transporting Ca2+-ATPase in the pellicles of Paramecium tetraurelia. Biochim Biophys Acta BBA—Biomembr. 1990;1029: 241–251. pmid:2147112
- 105. Husser MR, Hardt M, Blanchard M-P, Hentschel J, Klauke N, Plattner H. One-way calcium spill-over during signal transduction in Paramecium cells: from the cell cortex into cilia, but not in the reverse direction. Cell Calcium. 2004;36: 349–358. pmid:15451619
- 106. Sedar AW, Porter KR. The fine structure of cortical components of Paramecium multimicronucleatum. J Biophys Biochem Cytol. 1955;1: 583–604. pmid:13278368
- 107. Zhang P, Jana S, Giarra M, Vlachos PP, Jung S. Paramecia swimming in viscous flow. Eur Phys J Spec Top. 2015;224: 3199–3210.
- 108.
Graf B. Quaternions and dynamics. ArXiv08112889 Math-Ph. 2008 [cited 12 Jan 2021]. http://arxiv.org/abs/0811.2889
- 109. Moolenaar WH, De Goede J, Verveen AA. Membrane noise in Paramecium. Nature. 1976;260: 344–346. pmid:1256576