An electrophysiological and kinematic model of Paramecium , the 1 “swimming neuron” 2

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 17 triggered by a calcium-based action potential. For this reason, several authors have called Paramecium 18 the “ swimming neuron ” . Here we present an empirically constrained model of its action potential 19 based on electrophysiology experiments on live immobilized paramecia, together with simultaneous 20 measurement of ciliary beating using particle image velocimetry. Using these measurements and 21 additional behavioral measurements of free swimming, we extend the electrophysiological model by 22 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. 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.

. The avoiding reaction of Paramecium. A, Typical spontaneous avoiding reaction: the ciliate swims backward, then turns and eventually resumes forward swimming, while spinning around its main axis during top). Such stimuli are known to trigger different voltage-gated currents: a small inactivating calcium 126 current (Nakaoka and Iwatsuki, 1992;Preston et al., 1992aPreston et al., , 1992b, and a strong inward rectifier K + 127 current IKir (Oertel et al., 1978;Preston et al., 1990). We do not model the small calcium current, which 128 is related to the escape reaction, an increase in swimming speed triggered by hyperpolarizing stimuli.

delayed rectifier current measured in deciliated cells. A, Current-voltage relationship in
In the electrophysiological responses, we notice several differences with deciliated cells. First, an 241 upward deflection is apparent after stimulation, as illustrated in Fig. 4C (arrow). This deflection is due 242 to an inward current, the Ca 2+ current. This current can be estimated by subtracting the estimated leak 243 current from the capacitive current (passive properties estimated by model fitting, see below). With a 244 pulse of intensity I = 1.5 nA, we find that the inward part of that current peaks at about -2 nA (Fig. 4D).

245
This is an underestimation since part of the inward current may be masked by the K + current, but it is 246 comparable to previous estimations in voltage-clamp (Oertel et al., 1977). This current is known to 247 activate and inactivate quickly, within a few ms (Brehm et al., 1980;Oertel et al., 1977), as can be seen 248 on Fig. 4D. By integrating this current, we calculate that it should lead to an increase in intraciliary 249 calcium concentration of about ∫ /2 = 10 µM (a lower estimate, because of masking by K + currents),

250
where ≈ 1700 µm 3 is the estimated ciliary volume (see Methods, Electrophysiological modeling) and

251
F is the Faraday constant. This is well above the threshold for ciliary reversal, which has been 252 estimated at about 1 µM by exposing Triton-extracted cells to variable concentrations of calcium 253 (Naitoh and Kaneko, 1972).

254
In the electrophysiological response (Fig. 4A), we also observe small oscillations, due to the interplay 255 between Ca 2+ and K + currents, and a pronounced hyperpolarization after the pulse. This 256 hyperpolarization is due to a calcium-activated K + current IK(Ca). This current has been previously 257 characterized electrophysiologically (Saimi et al., 1983;Satow and Kung, 1980), as well as genetically

259
Thus, we included the following currents in our model: a leak current IL, a voltage-gated calcium 260 current ICa, with calcium-mediated inactivation, a delayed rectifier K + current IKd, and a calcium-261 activated K + current IK(Ca) (see Methods, Electrophysiological modeling). The calcium current ICa is 262 produced by ciliary channels similar to L-type Cav1.2 channels (Lodh et al., 2016). We modeled it 263 similarly to (Eckert and Chad, 1984;Standen and Stanfield, 1982), but we allow for several inactivation 264 binding sites (equations (14-16)). In addition, the current uses the Goldman-Hodgkin-Katz equation,

265
which is more appropriate than the linear driving force (ECa -V) when intra-and extracellular properties of plasma membrane calcium pumps (PMCA), also present in the cilia (Yano et al., 2015, Thus, calcium concentration is not directly measured, but indirectly constrained by several processes: voltage-clamp measurements (Oertel et al., 1977), the calcium current is transient and peaks at about 288 ~3.5 nA (Fig. 4E3). A residual current remains, so that calcium concentration remains high during 289 stimulation (Fig. 4E4). This is consistent with the fact that ciliary reversal can last for many seconds 290 when the membrane is depolarized (Machemer, 1974). With large currents, calcium concentration 291 raises to about 22 µM, similar to previous estimations.

292
The voltage-gated potassium current is delayed relative to the calcium current, and the calcium-293 activated K + current raises more slowly and is only dominant during repolarization (Fig. 4E3), with a 294 maximum of 0.65 nA. This is consistent with previous studies of that current (Satow and Kung, 1980). Figure 4F shows the three different currents during the action potential shown in Fig. 4F. As previously 296 argued, the calcium-activated K + current has a small contribution to the early current (Oertel et al., 297 1977).

300
Estimated conductance is not well constrained and often very large. This is presumably because the 301 peak current is mainly determined by the inactivation properties, and therefore the conductance 302 parameter is not well constrained.  Kung, 1980). There are about nKCa = 3 binding sites is (3.5 ± 2.3).

314
Finally, cilia revert at about 2.4 µM (log10(Kmotor in M) = -5.9 ± 0.8). This is close to measurements with 315 triton-permeabilized cells, reporting about 1 µM. We note that this and other concentration 316 parameters depend on the estimation of intraciliary volume, which is approximate.

317
Overall, parameters of the fitted models are compatible with known properties of the currents and of 318 ciliary reversal. 319 320 . CC-BY 4.0 International license perpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in observations of free swimming (Jennings, 1904;Machemer, 1972). During a pulse that triggers an 345 action potential, the flow is directed towards the anterior end, about 9° to the right (Fig. 5A, red). This

354
To understand the relationship between ciliary beating patterns and kinematic parameters, in 355 particular , we examine a spherical model of radius 60 µm (

360
At rest (low calcium concentration), cilia beat towards the rear, slightly to the right (Fig. 5E, first 361 column), so that the fluid produces a force towards the front, slightly to the left. If the direction of 362 ciliary beating is identical everywhere in spherical coordinates (that is, in terms of the cardinal 363 directions North/South/East/West, Fig. 9), then the force field over the sphere is symmetrical with 364 respect to the main axis. This makes both the total force and the total torque align with the main 365 (antero-posterior) axis, and therefore v and are also aligned with that axis, that is, = 0. The 366 organism then moves forward, with a spinning movement around the axis. We adjust the force so that 367 the velocity is 500 µm/s, which makes the sphere spin at about 1.8 Hz. Upon stimulation, when calcium 368 concentration is high, cilia revert and beat forward (Fig. 5E, second column), making the organism 369 move backward.

370
Thus, the organism cannot turn unless there is some asymmetry in the ciliary beating pattern.

382
sphere also spins about 4 times faster ( /(2 ) ≈ 7.2 Hz). It returns to moving forward when all cilia . CC-BY 4.0 International license perpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for this this version posted February 17, 2022. ; https://doi.org/10.1101/2022.02.15.480485 doi: bioRxiv preprint 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 387 the precise pattern of ciliary reversal empirically. Nevertheless, it is possible to demonstrate that cilia 388 revert asynchronously. We measured particle flow separately in two regions of the field, beyond the quality was sufficient (see Methods, Particle image velocimetry). We measured the mean angle of the 391 flow field during stimulation with current pulses between 1 nA and 5 nA; weak pulses were not 392 included because particle density tended to be lower due to sedimentation (weak pulses were 393 recorded after strong pulses).

424
. CC-BY 4.0 International license perpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in

459
The characteristics of the avoiding reaction also depend on stimulus duration (Fig. 6G). If the pulse 460 amplitude is fixed (I = 0.1 nA) and its duration is increased, then the duration of backward swimming

463
When we examine spontaneous avoiding reactions of freely swimming paramecia, we find that both 464 backward swimming duration and reorientation angle vary broadly (156 ± 81 ms and 114 ± 66 °, 465 respectively) (Fig. 6H), and there is a small although highly significant correlation (linear regression,

472
. CC-BY 4.0 International license perpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in

481
First, we consider an organism swimming towards a generic object, which triggers a depolarizing 482 current when in contact with the membrane (for example a chemical substance, or hot water) (Fig. 7).

483
Thus, we simply consider that the stimulus current is proportional to the surface area in contact with 484 the stimulus (see Methods, Sensory transduction). In contrast with previous situations, the stimulus is 485 not pre-determined but depends on behavior. As Dewey pointed out (1896), "the motor response 486 determines the stimulus, just as truly as sensory stimulus determines the movement.". When the organism 487 touches the object, a current is triggered, which depolarizes the membrane (Fig. 7A, 7B). As the cell 488 swims into the object, the current increases until an action potential is triggered. The cell then swims 489 backward, moves out from the object and the current stops. Thus, the sensory current is necessarily 490 small and short, because the organism withdraws as soon as the current reaches threshold. This results 491 in a small avoiding reaction, and the organism bumps again repetitively against the object until it 492 finally escapes (Movie 1).

493
Larger movements can be obtained if sensory transduction has slower kinetics (Figs. 7C, D and Movie 494 2). Here, the sensory current follows the stimulation with a time constant of 40 ms, modelled with first 495 order kinetics (to simplify, the spatial spread of channels is not modeled; equation (24)). In this case,

496
. CC-BY 4.0 International license perpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for this this version posted February 17, 2022. ; https://doi.org/10.1101/2022.02.15.480485 doi: bioRxiv preprint 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. the organism spins while swimming, so that its oral groove takes a different position when it touches 501 the object.

502
We now consider that the object is a disk in a square environment ( Fig. 8A and Movie 4). To avoid 503 boundary effects, we consider that the environment has the topology of a torus (paramecia escaping 504 to the left reappear to the right). To account for spontaneous avoiding reactions (occurring at a rate of 505 about 0.2 Hz in our behavioral measurements), we added a noisy current to the membrane equation.

506
The paramecia are modeled as in Figs   be seen that paramecia tend to aggregate inside the disk, mostly near the boundary. This is due to the 525 small avoiding reactions but also to the curvature of the disk.

526
. CC-BY 4.0 International license perpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in  (Berg, 1975).

561
The electrophysiological model was built by model fitting to current clamp data, with calcium-562 dependent properties indirectly constrained by ciliary reversal data. This method recovered 563 properties of individual currents compatible with previous measurements obtained by different 564 means. For example, the calcium threshold for ciliary reversal was estimated to be ~2 µM, the same 565 order of magnitude as measured by varying extracellular calcium concentration in Paramecium with 566 permeabilized membranes (Naitoh and Kaneko, 1972). This is notable since calcium was not measured 567 but only inferred from electrophysiology (indirectly through the estimation of the calcium current by 568 the fitting procedure). The fitting procedure also determined that the calcium-dependent K + current is 569 small during the action potential, as previously determined with voltage-clamp experiments (Oertel et

581
The integrated model shows helicoidal swimming with graded avoiding reactions, where backward 582 duration swimming and reorientation angle increase with stimulus strength or duration. As observed 583 in spontaneous behavior, the model can also slightly reorient without swimming backward.

584
Behavior of the autonomous model is more complex than stimulus-response experiments, because the 585 relation between sensory stimulus and motor response is circular. In particular, we noticed that the 586 interaction with an object (e.g. a chemical substance) critically depends on the properties of sensory

614
Another limitation is we did not measure ciliary beating directly, but rather its effect on the fluid. This

619
. CC-BY 4.0 International license perpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in

654
Physiologically, asynchronous ciliary reversal can be due to differential calcium sensitivity, that is, the 655 calcium threshold for reversal might vary across cilia. This is the implicit assumption of our model. It  Kung, 1985). The model included two calcium channels (fast and slow), for which there is no 670 electrophysiological support, and their relative activation was given as a function of time after stimulus 671 start (i.e., it is not modeled). The calcium-dependent K + current was not included. Neither 672 electrophysiological model was fitted to data, and neither was coupled to a kinematic model.    CC-BY 4.0 International license perpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for this this version posted February 17, 2022. ; https://doi.org/10.1101/2022.02.15.480485 doi: bioRxiv preprint al., 1986). The tube was shaken before collecting cells to perform an experiment.

705
The culture method differed slightly for the behavioral measurements with freely swimming an experiment, bacteria were first grown in 5 mL of WGP for 24 h at 27°C, then paramecia were grown  (Iwatsuki andNaitoh, 1983, 1982). Movies of the swimming paramecia are 200 s long (see Movie 7).

734
To limit hard drive space, images are stored without their background with lossless compression (TIFF

745
. CC-BY 4.0 International license perpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in

815
In each frame, we calculated the mean angle of the velocity vector over the entire field, using circular 816 mean (argument of the mean complex unit vector; occasional missed frames were discarded). We then 817 subtracted the angle of the anteroposterior axis. The position of anterior and posterior ends was 818 measured manually. As the two ends can be visually ambiguous, they were automatically corrected (by 819 swapping) when the flow measured before stimulation was directed towards the anterior end 820 (indicating backward swimming).

821
In Figure 5A, for each cell we averaged the mean angle over all currents and over the 300 ms before 822 stimulus (blue) or over the second half of the stimulus (red), for positive currents (<5 nA).

823
We also calculated the mean angle in the anterior and posterior regions as indicated in Fig. 5F. Each 824 region is a half-plane orthogonal to the main axis, starting at one end. For this analysis, we selected n 825 . CC-BY 4.0 International license perpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in

855
where C is membrane capacitance, : = : ( : − ) is the leak current, IKd is the delayed rectifier K + 856 current responsible for repolarization, ICa is the ciliary voltage-dependent Ca 2+ current, IK(Ca) is the 857 calcium-activated K + current, IKir is the inward rectifier K + current and I is a stimulating current. We 858 did not include a few other electrophysiologically identified currents that are less relevant for this 859 study, namely: Na + (Saimi, 1986;Saimi and Ling, 1990) and Mg 2+ (Preston, 1998, 1990) currents (since 860 . CC-BY 4.0 International license perpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for this this version posted February 17, 2022. ; https://doi.org/10.1101/2022.02.15.480485 doi: bioRxiv preprint calcium currents responsible for the escape reaction (Nakaoka and Iwatsuki, 1992; Preston et al., where p = 1 or 2 (p = 2 in the final version). We made this simple modeling choice because this current 868 was only used as a way to infer the reversal potential EK. In particular, we did not include inactivation

872
For the delayed rectifier current IKd, we tested several models of the type: We tested two classes of models. The Hodgkin-Huxley model is:

884
The Boltzmann model is:

886
. CC-BY 4.0 International license perpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for this this version posted February 17, 2022. ; https://doi.org/10.1101/2022.02.15.480485 doi: bioRxiv preprint

887
The voltage-gated calcium current ICa is a calcium-inactivated current located in the cilia (Brehm et al.,

908
This is similar to the model of (Standen and Stanfield, 1982), except that the number of sites nCa is 909 allowed to be greater than 1, because we found that this was necessary to fit our data.

910
A calcium-activated K + current has been identified by comparison with Pawn mutants lacking 911 functional voltage-activated calcium currents (Satow and Kung, 1980

917
. CC-BY 4.0 International license perpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for this this version posted February 17, 2022. ; https://doi.org/10.1101/2022.02.15.480485 doi: bioRxiv preprint (Iwadate, 2003;Klauke and Plattner, 1997). We chose [Ca 2+ ]0 = 100 nM. Calcium enters the cilia when calcium pumps (PMCA) have been identified in the basal membrane with low affinity, around 10 -7 M 926 (Wright and van Houten, 1990), and also in the cilia (Yano et al., 2015(Yano et al., , 2013. Suppressing the ciliary 927 PCMAs by RNA interference prolongs backward swimming, which means that they are involved in the 928 removal of calcium after an action potential. In principle, calcium can also diffuse to the basal cytosol.

929
However, this has not been observed (Husser et al., 2004). This might be because of cilia volume 930 compared to the cell, or because calcium is buffered at the base of cilia. Both phenomena can be 931 modeled by diffusion to the cilium base, with fixed resting concentration at the boundary.

932
We lump these diverse mechanisms into two simple processes: a linear process, with rate proportional

937
where v is the volume of cilia and F is the Faraday constant. It can be seen that the role of the high-

938
affinity process in this model is to counteract the calcium flow at rest, namely = #9MN / cilia , while 939 the low-affinity process independently tunes the rate of calcium removal after an action potential.

944
We used the previous estimate cilia = 1700 µm 3 from (Oertel et al., 1977), which is compatible with 945 these bounds, but the uncertainty is large. In addition, the effective volume might be smaller because

955
. CC-BY 4.0 International license perpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in

974
We set min = 13° based on our measurements, and max = 90°, to allow for planar rotations as  Table 1.  Fig. 3, we fitted the models described above with IL and IKd, with least square minimization of the 1041 error on the reading electrode potential, taken from pulse start to 50 ms after the pulse. The stimuli 1042 were 100 ms pulses with amplitude between 0 and 4 nA in 300 pA increments. Parameters EK and C 1043 were taken from the previous fit to hyperpolarized responses. To make sure the short onset is well each window is equally weighted (meaning that a data point in the first window contributes more than a data point in the second window).

Mean
The local force is tangent to the sphere and oriented with an angle , where = 0 is the direction of 1112 the meridian, pointing South. Cartesian coordinates are chosen so that the x axis is dorso-ventral, the

1117
. CC-BY 4.0 International license perpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The local torque is = × where is a radius. In Cartesian coordinates, we obtain:

1128
Three different patterns are represented in Fig. 5E. In the forward pattern, the local angle is uniform: 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 N^N and N^N are aligned with the z axis (anteroposterior).

1133
This can be seen in the formulas above by integrating with respect to , which yields 0 for the x and y 1134 coordinates.

1135
In  1146 . CC-BY 4.0 International license perpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in We consider that the organism is an object moving by rigid motion. The organism is characterized by reference frame is chosen as in the spherical model above, so that z>0 points towards the anterior end  1177 1178 . CC-BY 4.0 International license perpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in In all simulations, the electrophysiological model fitted to the cell shown in Fig. 4

is integrated with
In Fig. 6D, models are simulated in a plane with 2 ms pulse currents of amplitude 0.3, 0.5 and 5 nA, at 6G, the stimulus is a 100 pA pulse of duration 0 to 100 ms. The reorientation angle is calculated as the 1188 change in angle before and after the stimulus, averaged over 1 second (the spinning period).

1190
Interaction with a stimulus 1191 In Fig. 7 and all subsequent figures, trajectories are constrained to a plane. In Fig. 7A and B, the stimulus 1192 is a half-plane. It produces an instantaneous depolarizing current, proportional to the fraction of the 1193 cell shape inside the stimulus (see Sensory transduction), with maximum 5 nA. In Fig. 7C and D, the 1194 stimulus additionally goes through a low-pass filter with time constant 40 ms, representing the activation/deactivation rate of the receptors.
In Fig. 8A and 8B, the stimulus is a disc of radius 1 mm within a 4 mm torus.  In Fig. 8C, the environment has the topology of a cylindrical surface, i.e., circular in the small dimension 1211 (500 µm) and linear in the long dimension (25 mm). The stimulus is a linear gradient of 100 pA/mm, 1212 with transduction modeled with adaptation as for the attracting disc (Fig. 8B).

1213
1214 . CC-BY 4.0 International license perpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in In Fig. 8D, cells produce CO2 by breathing, which then diffuses and acidifies the fluid. This is modelled coefficient of CO2. In water at 25°C, ≈ 0.002 mm 2 /s but we accelerate it by a factor 5 to speed up the Brette R. 2021. Integrative Neuroscience of Paramecium, a "Swimming Neuron." eNeuro 8:ENEURO.0018-

1248
. CC-BY 4.0 International license perpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in

1286
. CC-BY 4.0 International license perpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in

1325
. CC-BY 4.0 International license perpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in