Chemophoresis engine: A general mechanism of ATPase-driven cargo transport

Cell polarity regulates the orientation of the cytoskeleton members that directs intracellular transport for cargo-like organelles, using chemical gradients sustained by ATP or GTP hydrolysis. However, how cargo transports are directly mediated by chemical gradients remains unknown. We previously proposed a physical mechanism that enables directed movement of cargos, referred to as chemophoresis. According to the mechanism, a cargo with reaction sites is subjected to a chemophoresis force in the direction of the increased concentration. Based on this, we introduce an extended model, the chemophoresis engine, as a general mechanism of cargo motion, which transforms chemical free energy into directed motion through the catalytic ATP hydrolysis. We applied the engine to plasmid motion in a ParABS system to demonstrate the self-organization system for directed plasmid movement and pattern dynamics of ParA-ATP concentration, thereby explaining plasmid equi-positioning and pole-to-pole oscillation observed in bacterial cells and in vitro experiments. We mathematically show the existence and stability of the plasmid-surfing pattern, which allows the cargo-directed motion through the symmetry-breaking transition of the ParA-ATP spatiotemporal pattern. We also quantitatively demonstrate that the chemophoresis engine can work even under in vivo conditions. Finally, we discuss the chemophoresis engine as one of the general mechanisms of hydrolysis-driven intracellular transport.


Introduction
Cell polarity regulates the direction of intracellular transport for cargos, such as organelles and macromolecules, by taking advantage of chemical gradients sustained with the aid of ATP or GTP hydrolysis [1].For example, it is well known that eukaryotic cell polarity factors, such as Rho, GTPase, and Cdc42, regulate the orientation of cytoskeleton members so that molecular motors can carry cargo directionally on the cytoskeleton, contributing to cell movement [2], cell growth [3], and axon guidance [4].Although the transport by the cytoskeleton is one of the most commonly observed mechanisms, the transport directly mediated by chemical gradient, if its existence is confirmed, should be of importance as a general mechanism for cargo transport as well, which we refer to as cargo chemotaxis here.
Similar to the Min system, the pole-to-pole oscillation of ParA [14,23,[63][64][65][66][67] also emerged in the ParABS system through the stimulation of ParA ATPase activity by ParB on the PC [13][14][15][22][23][24].Interestingly, a plasmid chases a ParA focus, following its oscillatory movement along the long host-cell axis [23], leading to oscillatory motion.In a recent in vitro experiment mimicking a ParABS system, Vecchiarelli et al. elegantly demonstrated the formation of directed motion of a cargo corresponding to a plasmid, referred to as "cargo surfing on ParA-ATP traveling wave" [28][29][30].Hence, for both Min and Par systems, the emergence of traveling waves and the pole-to-pole oscillation of ATPase have been reported.The mechanism driving the plasmid motion, however, remains elusive [18], whereas the pattern dynamics of the Min system can be described by well-defined reaction-diffusion equations [58][59][60][68][69][70][71][72].
In the present study, we extend our previous model and propose a chemophoresis engine as a general mechanism of cargo motion, which transforms chemical energy into directed motion via self-organization of the traveling wave, and then apply it to the plasmid motion in a ParABS system.In the previous study, the plasmid was assumed to be a point particle, where static equi-positioning and symmetric ParA-ATP distribution were robustly maintained [43].However, such model with a zero-size limit is unrealistic, considering intracellular dynamics [26,73] or reconstructing in vitro experiments performed by [28].Here, by considering the finite size of plasmids explicitly, we show that organization of directed motion is possible via spontaneous symmetry breaking in the ParA-ATP pattern.We then recapitulate plasmid positioning to better describe the spatiotemporal profiles of ParA-ATP concentration and movement of plasmids.Actually, in the model presented here, the net chemophoresis force acts on the plasmid (PC) through the concentration difference between its ends, which is self-sustained by the high ATP hydrolysis stimulation.This self-driven mechanism leads to the directed motion of plasmids, as well as their equi-positioning [16,24,64,67], and pole-to-pole oscillation as observed in bacterial cells and in vitro experiments [23,[64][65][66][67].We mathematically show the existence and stability of the plasmid-surfing pattern, which allows cargo-directed motion through the symmetry-breaking transition of the ParA-ATP spatiotemporal pattern.We also indicate that plasmid size is a relevant parameter for the emergence of its directed movement.Finally, we quantitatively validate that the chemophoresis engine can work with parameters capturing in vivo conditions.

Chemophoresis force
First, we briefly reviewed the chemophoresis force, a thermodynamic force acting on the cargo in the direction of the increased concentration of a chemical that can be bound on the cargo (See S1 Text for details.)We considered that a cargo was placed and moving in a d-dimensional space r 2 R d .The cargo had N molecular sites B, on each of which m molecules of chemical X was bound to form a complex Y at position r = ξ.At each site, the reaction mX(ξ) + B ⇄ Y occurred and was at chemical equilibrium.If a spatial gradient of chemical concentration X exists, the cargo is thermodynamically driven in the direction of the decreased free energy or the increased chemical potential of X [42,43].Here, such a gradient of the chemical potential μ(r) was assumed to be sustained externally through several active processes, supported by spatially distributed chemical gradients.We referred to the phenomenon as chemophoresis.The formula of the chemophoresis force was: Here, mðrÞ ¼ � m þ k B T ln xðrÞ, where x(r) is the concentration of X, K d is the dissociation constant, and m is the number of binding molecules corresponding to the Hill coefficient of the reaction.With this chemophoresis force, the cargo moved in the direction such that the concentration x(r) increased even under thermal fluctuation ( [42,43], S1 Text).For the force to work, the reaction mX + B ⇄ Y was required to reach chemical equilibrium fast enough for cargo motion.Therefore, we showed that the chemophoresis force was one of the fundamental thermodynamic forces driven by physicochemical fields.Note that the force had an entropic origin from the viewpoint of statistical mechanics.See also Ref. [43] for details of the derivation from the viewpoint of thermodynamics and statistical mechanics.
To understand the origin of chemophoresis, it should be noted that microscopic binding events of X do not directly generate the force.Rather, the force works in the direction of larger frequency of the binding events (or larger time fraction of binding states) that was realized in the spatial location with a larger concentration of molecules in a chemical bath.The chemical gradient biases the binding frequency of X in a space-dependent manner.In other words, chemophoresis is driven by general thermodynamic force as a result of the free-energy (entropy) difference.It can also be derived by coarse-graining microscopic processes (S1 Text), whereas the macroscopic derivation implies its generality independent of specific microscopic models [35][36][37].On the other hand, both the macroscopic (thermodynamic) and microscopic (statistical physics) theories are equivalent to each other, in that the force is generated with the aid of spatial asymmetry of molecule numbers bound on the bead, if its radius is finite.Further, for the chemophoresis force to act, X molecules do not necessarily have to bind cooperatively to the bead (as in the case of m = 1); if the concentration gradient of bound molecules is generated, the resultant free energy difference between its ends leads to the net chemophoresis force.

Chemophoresis engine for plasmid partition
We then applied the chemophoresis formula to plasmid motion.As the reaction on the cargo consumed chemical X, its concentration changed; therefore, studied its RD equation.It was introduced for ParA-ATP ( [42,43], S1 Text), which recapitulates the central-and equi-positioning of plasmids [42,43].It was also adopted successfully to explain the directed movement of beads in an in vitro experiment by Vecchiarelli et al. [28].We considered a plasmid i(1 � i � M) placed into and moving in a d-dimensional space r 2 R d (d = 1 or 2) (Fig 1A).ParA-ATP dimers were bound to a PC on plasmid i at position r = ξ i .m ParA-ATP dimer molecules interacted with ParB, which stimulated ParA ATPase activity at a catalytic rate k [74]; N ParB molecules were assumed to be recruited to each PC at r = ξ i .Because ParA could not bind PC when it was not combined with ATP, free ParA products were released from the PC immediately after ATP hydrolysis.Thus the reaction was presented as follows: Through this reaction on PC i at r = ξ i , each plasmid acted as a sink for ParA-ATP and induced a concentration gradient of this protein.In the early model, the size of plasmids was assumed to be zero ( [42,43], S1 Text).However, the model is still too unphysical to better reconstruct the movement of plasmids with a finite size in bacterial cells [26,73] as well as that of micro-sized beads in in vitro experiments [28].To better describe spatiotemporal profiles of ParA-ATP concentration and directed movement of the plasmids, we considered plasmids (or PCs) as spheres with a radius of l b whose value is reported to be l b � 0.075 μm in bacterial cells according to [26,73] and l b = 1.0 μm in in vitro experiments [28].Here, in order to discuss general situations, the derived equations were first rescaled by a dimensionless form and then numerical simulation was performed.Denoting the dimensionless concentration of ParA-ATP dimers on a nucleoid as u(r), the normalized RD equation was written as follows (see S1 Text for its derivation): where the first and second terms represent the diffusion of ParA-ATP and its chemical exchange at a normalized constant rate with the cytoplasmic reservoir (denoted by its normalized concentration), respectively.The last term denotes the inhibition by ParB on the M PCs.K d is the normalized dissociation constant of the reaction mX + B ⇄ Y, and m is the Hill coefficient.V is the d-dimensional volume of the bead with a radius of l b .χ = kN/V is a maximum rate for ParA-ATP hydrolysis by ParB on each PC (S1 Text).Furthermore, θ(r) is a step function representing the space each PC occupies to describe the hydrolysis reaction space.Only within |r − ξ i | < l b , the reaction occurred.Without the last term (if χ = 0), u(r) reached a homogenous equilibrium state, u (r) = 1.In contrast, the normalized equations of motion for plasmids were represented as follows: with thermal noise hη i (t)i = 0 and hη i ðtÞ � η j ðt 0 Þi ¼ 2dDd ij dðt À t 0 Þ, and ε ¼ DN=V.Here,

Chemophoresis engine captures observed plasmid dynamics
Chemophoresis engine can recapitulate equi-positioning, directed movement, and poleto-pole oscillation.First, we considered the motion in a one-dimensional (1D) space (d = 1), that is, on a nucleoid matrix along the long cell axis where a plasmid i(1 � i � M) was positioned at x = ξ i 2 [0, L] (Fig 1B ); The Neumann boundary condition was adopted for the RD equation: ru(0) = ru(L) = 0. To confine the plasmids to the host cell x 2 [0, L], we placed the reflection walls at x = 0 and x = L.This could be explicitly represented as In general, the plasmid showed directed motion for a larger χ (= maximum rate of ParA-ATP hydrolysis).This result was plausible because the larger χ generates the sharper gradient of ParA-ATP, leading to the larger chemophoresis force to enable the persistent directed motion of the plasmid.
Similarly, for M > 1 cases, plasmid dynamics qualitatively changed among stochastic switching, steady equi-positioning, and directed movement followed by an oscillatory mode as χ increased (S1 and S2 Figs).We then examined how plasmid dynamics switched from a static to an oscillatory mode with increasing χ.The switch in plasmid dynamics occurred through a symmetry-breaking transition of the ParA-ATP spatiotemporal pattern.Interestingly, the oscillatory behavior of plasmids did not disrupt time-averaged equi-positioning.Steady multimodal distribution of plasmids was sustained (S1 and S2 Figs).As reported previously [42,43], the regular positioning of plasmids is due to the effective inter-plasmid repulsive interaction derived from the chemophoresis force.The plasmid acting as a sink for ParA-ATP contributed to the formation of a concentration gradient, which increased with the distance from the plasmid.Other plasmids were subjected to the chemophoresis force caused by the gradient in the direction of increasing ParA-ATP concentration so that they were forced away from the former.The former plasmid was also subjected to the chemophoresis force caused by the gradient derived from the latter, resulting in mutual repulsion among the plasmids.The mutual repulsive interaction contributed to the robustness of the positional information generated by the chemophoresis engine.This plasmid separation scenario by such repulsive interactions is consistent with a previous observation [16].
Chemophoresis engine mathematically validates plasmid surfing on the traveling ParA-ATP wave.In a recent report, [28], Vecchiarelli et al. demonstrated the directed movement of micro-sized beads that mimic plasmids.A theoretical explanation using more realistic model with finit-sized plasmids is needed, as the previous models [43] assumed vanishing-size plasmids; Hence, we introduced the l b -sized plasmids in Eqs 3 and 4, to discuss a possible symmetry-breaking transition more realistically.Here, we analytically examined plasmid surfing on the ParA-ATP traveling wave, focusing only on M = 1, without thermal fluctuation, for a 1D case under a periodic boundary condition.
S3 Furthermore, we performed linear stability analysis of surfing-on-wave solution for S3 and S4 Eqs in S2 Text against external disturbances, in order to examine if tiny perturbation Plasmid size is a relevant parameter for the emergence of its directed movement.In the limit with vanishing plasmid size l b !0, the mathematical form of the present model is reduced to our earlier model (S1 Text).However, whether any spontaneous directed motion emerges in the limiting case or not was not discussed in our previous report [43].In such vanishing size limit, χ and ε diverge: χ = kN/(2l b ) ! 1 and ε ¼ DN=ð2l b Þ !1; nevertheless steady surfing-on-wave solutions can exist in a certain range of k; D; and N. By using the analytical solution obtained above (S7-S9 Eqs in S2 Text) and by taking a large system size limit L ! 1, the steady velocity in the limit is analytically obtained as (S25 Eq in S2 Text): v ¼ � ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi black dot plots).The numerical results showed slight deviation from the analytical solution for the ranges of small and large l b , suggesting the breakdown of the approximation u(x) � K d .In general, the directed motion with larger l b requires a larger k and N.This is because a plasmid (or bead) with large l b cannot generate sufficient concentration gradient and the resultant chemophoresis force is not sufficient to realize its self-driven directed motion.Finally, these results indicate that plasmid size is a relevant parameter for the emergence of its directed movement.Actually, in a previous mathematical study [37], its size relative to cell length was reported to play a different role in plasmid partition than the present result.

Chemophoresis engine can work in in vivo conditions
From the normalized equations for chemphoresis engine (Eqs 3 and 4), we here confirmed that plasmid surfing on the traveling wave emerged through a symmetry-breaking transition with the change of χ and ε (Fig 2).However, it remains unclear whether the above results can be quantitatively reproduced for the parameter values reported experimentally in bacterial cells.To quantitatively confirm that the directed motion is self-organized by the chemophoresis engine, we ran the simulations again with parameters capturing in vivo conditions (S1 Table ), as described below.
A recent study suggests that ParB localizes around plasmids by a ParA-ATP dependent phase separation mechanism to form PCs as droplets, accompanied by an enhanced hydrolysis activity of ParA-ATP through the increase in local density of ParB within the PCs, finally leading to a successful plasmid partition [75,76].As a simple demonstration, we carried out the simulation by simultaneously changing the maximum rate of ParA-ATP hydrolysis χ = kN/ (2l b ) and the strength of the chemophoresis force ε ¼ DN=ð2l b Þ under the control of the local density of ParB on a PC, ρ ≔ N/(2l b ).We increased the cell length, which may reflect the cell cycle progression, whereas we adopted the parameters capturing in vivo conditions listed in S1 Table .We changed ρ and the cell length = L × l with L as a normalized system size, scaled by l ¼ ffi ffi ffi ffi ffi ffi ffi ffi demonstrating that the chemophoresis engine can recapitulate plasmid surfing even in in vivo conditions.Interestingly, the bifurcation points differ by different cell lengths; that is, the emergence of directed movement can be cell-length dependent.This bifurcation to the directed motion against the cell length occurs for both M = 1 and M = 2 (S7B1 and S7B2 Fig) .In addition, the bifurcation diagram for M = 1 coincides with that for M = 2 when plotted as a function of (cell-length)/M for the both cases (S7B Fig) , suggesting that controlling the inter-plasmid distance is important for the emergence of directed motion.Such a scaling for bifurcation diagram among different values of M, as was mentioned in [37], can enable robust equi-positioning of plasmids during cell elongation [16,67].

Discussion
In this study, to consider a generalized model of the plasmid partition ParABS system, a chemophoresis engine was introduced as a coupled dynamical system among the equations of motion for plasmids and the RD equation for ParA-ATP (Fig 1A).In the model, plasmid dynamics switched from static to dynamic mode with an increase in the maximum rate of ATP hydrolysis χ.The engine demonstrated equi-positioning, directed movement, and poleto-pole oscillation, as observed in bacterial cells and in vitro experiments (Figs 1C, S1 and S2).Note that despite the plasmids' oscillatory behavior, the regular positioning distributions were sustained (S1 and S2 Figs) due to an effective inter-plasmid repulsive interaction derived from the chemophoresis force, indicating the robustness of positional information generated by the chemophoresis engine.By simplifying Eqs 3 and 4, and introducing a space-time coordinate, we mathematically showed the existence ( Although we analyzed the existence and stability of the surfing-on-wave pattern only in a noiseless situation (Fig 2 ), plasmids (or cargos) are always subjected to thermal fluctuations in cellular environments.Then, the plasmid motion was described by Langevin equation Eq 4. Further, we needed to elucidate that the plasmid-surfing-on-traveling-wave pattern remains robust against thermal fluctuations.For the chemophoresis force to act effectively, the force must be larger than the thermal noise, as discussed in a previous report [43].Any force weaker than thermal noise cannot sustain even regular positioning [43].We also examined how equipositioning of plasmids can overcome thermal noise disturbances in a 1D case (S1 and S2 Figs).We confirmed that plasmid location dynamics shows a transition from stochastic switching to (freezing) steady equi-positioning (S1 and S2 Figs) as χ is increased, finally leading to persistent directed motion of the plasmids.This result suggested that the chemophoresis force dominates and directed movement of plasmids can overcome against thermal fluctuation for large hydrolysis rate.
We propose a chemophoresis engine, a general mechano-chemical apparatus driving the self-motion of the intracellular cargo, as a means to elaborate the physical principles of ATPase-driven cargo transport [48][49][50][51][52][53].The engine is based on 1) a chemophoresis force that allows motion along an increasing ATPase(-ATP) concentration and 2) an enhanced catalytic ATPase hydrolysis at the cargo positions.ATPase-ATP molecules are used as fuel to supply free energy by applying the chemophoresis force along the concentration gradient, whereas cargos generate a concentration gradient by catalyzing the hydrolysis reaction on their surface.Note that each cargo, as a catalyst, does not consume ATP, but only modulates the concentration pattern.Through the coupling and synergy between 1) and 2), directed movement of the cargo is self-organized, showing a "surfing-on-traveling-wave" pattern (Fig 2C).The chemophoresis engine is based only on these two general mechanisms and is expected to explain how the transportation of diverse cargos in bacterial and eukaryotic cells is organized.Although we have focused on the gradient generated by the regulation of ATPase, the regulation of the concentration gradient via phosphorylation-dephosphorylation reactions is ubiquitous.Therefore, the chemophoresis engine resulting from the regulation of the hydrolysis of other factors, such as GTPase, should work for a variety of intracellular processes [77][78][79][80].
Our theory is derived from macroscopic thermodynamics under nonequilibrium conditions, and although we have applied it here to the ParABS system, it is general enough to be independent of individual microscopic models constructed for each molecular mechanism.In the present study, the mathematical model of the chemophoresis engine is constructed only by extracting the essential parts of the phenomena, so it can be applied directly to other systems with common reaction mechanisms such as hydrolysis.Indeed, it has been reported in in vitro experiments that the directed motion of a micro-sized bead is self-driven by the RNA gradient which RNA hydrolysis on the bead generates [81]; The authors later termed this phenomenon "autochemophoresis" [82].
Furthermore, recent studies demonstrated substrate-driven chemotactic behaviors of metabolic enzymes (single-molecule chemotaxis) [83,84], and discussed a mathematical model to recapitulate experimental results in the subsequent study [85].Interestingly, the authors proposed the exact same thermodynamic mechanism as the chemophoresis force described previously [42,43].Therefore, we expect that the chemophoresis engine can also be applied to self-chemotactic behaviors even at a single-molecule level even though in the present study, self-chemotaxis is applied to a cargo size ranging from 50 nm to 1μm.However, to describe nanoscopic chemotaxis, we need to extend thermodynamics of chemophoresis to a stochastic one which is valid even under large thermal/chemical fluctuations.
The merits of the chemophoresis engine are as follows: Self-generated chemical gradient for the chemophoresis force to apply; not requiring a large space to maintain the external chemical gradient.The chemophoresis engine can be effective in a moderate space.Therefore, the chemophoresis engine would work for eukaryotic intra-nuclear processes by restricting the mobility of chemicals on a nuclear membrane or a nuclear matrix functioning as a scaffold matrix.Based on the generality of the chemophoresis engine as well as suggestive reports in other systems [81,82,86], we can apply the mechanism to other hydrolysis events, RNAs, receptors, and others.We propose the chemophoresis engine as a general mechanism for hydrolysisdriven cargo transports in cells.

Numerical methods for solving evolutionary equation, self-consistent equation, and eigenvalue equation
Evolutionary equations, Eqs 3 and 4 were computationally solved as a hybrid simulation between reaction-diffusion equation and Langevin equation.Euler scheme for Eq 3 and Euler-Maruyama scheme for Eq 4 were used as numerical algorithms.Real-valued self-consistent equation, S9 Eq in S2 Text and complex-valued eigenvalue equation S16 Eq in S2 Text were solved by using Newton-Raphson method.
Since ParA-ATP always exists as a dimer on a nucleoid, it is reasonable to simply consider a hydrolysis reaction without any cooperativity in a spatially limited space around a plasmid.Therefore, the Hill coefficient m was fixed as m = 1 in all the simulations.A normalized radius of plasmid l b was assigned to l b = 0.2 through the simulation except for Fig 2D.From tiny perturbations around the stationary solution U(z, t) = U st (z) + e λt δU(z), and z ξ (t) = e λt δz ξ , the eigenvalues λ were computed as a function of χ.The stability of the plasmid-surfing pattern (blue line) and the instability of the plasmid-localized solution (red line) was confirmed for χ > χ c in the absence of any positive real parameters of λ, Re[λ (v)] > 0, v 6 ¼ 0. It seems to regain stability at χ larger than χ � 5.9.However, such localized solution cannot be numerically realized for χ > 5.5 because U st (0) becomes negative and the solution is unphysical at χ � 5.5 (inset) as a result of breaking the approximation u(x) � K d .Therefore, the stable localized solution does not exist for w ≳ 5:5. is symmetrical, and the plasmid (red circle) is located at its minimum.In contrast, for the latter (B), the symmetry of u(r, t) is broken, indicative of a traveling wave.The minimum of the asymmetrical u(r, t) is positioned at a location shifted from where the plasmid lies, suggesting that the plasmid in "surfing" on the traveling wave."K d = 0.001, l b = 0.2, ε = 5, and L 2 = 10 × 10. (PDF) S1 Table .Model parameters capturing in vivo conditions.(PDF)

Fig 1 .
Fig 1. Chemophoresis engine can recapitulate equi-positioning, directed movement, and pole-to-pole oscillation.(A) Schematic representation of the chemophoresis engine.A plasmid moves in a d-dimensional space r 2 R d (d = 1 or 2).ParA-ATP dimer (green sphere) binds a partition complex (PC, magenta sphere) on the plasmid at position r = ξ i .ParA-ATP dimer molecules interact with ParB molecules (white spheres), which stimulate ParA ATPase activity at a catalytic rate.Because ParA cannot bind PC when it is not combined with ATP, free ParA products (blue sphere) are released from the PC immediately after ATP hydrolysis.Through this reaction on PC i at r = ξ i , each plasmid acts as a sink for ParA-ATP and induces a concentration gradient of this protein.(B) One-dimensional case, on a nucleoid matrix along the long cell axis where a plasmid i(1 � i � M) is positioned at x = ξ i 2 [0, L]. (C) The dynamics change among thermal motion, steady center-positioning, and directed movement followed by oscillatory mode as χ increases among χ = 0.5 (C1), χ = 2.5 (C2), and χ = 10 (C3) (two inner figures).(C1) The plasmid slightly tends to be localized at the cell center but it is still dominated by thermal fluctuations for M = 1 and χ ≔ kN/V = 0.5.(C2) It is stably localized at the cell center for M = 1 and χ = 2.5, and (C3) it shows directed movement, reflection at the end walls, and pole-to-pole oscillation for M = 1 and χ = 10.The corresponding ParA-ATP pattern dynamics also change among stochastic, steady center-positioning, and oscillatory waves (left).The oscillatory behavior of plasmids does not disrupt time-averaged center-positioning, but steady center-positioning of plasmids are sustained (Compare (C2) right and (C3), right).K d = 0.1, ε = 5, and L = 5.The distributions (right) were generated using 10 7 samples over 10 5 time step.https://doi.org/10.1371/journal.pcbi.1010324.g001 the relative diffusion coefficient of the plasmid to that of ParA-ATP (see S1 Text for details).The parameters to be assigned to Eqs 3 and 4 are K d ; l b ; m; M; N; D; k, and the system size L(= cell length).
40 through the simulation, and then we examined how the plasmid dynamics change with χ ≔ kN/V, where χ is the normalized maximum rate for ParA-ATP hydrolysis by ParB on each PC.For M = 1, the dynamics of u(x) and the plasmid, as well as the distribution of the plasmid position, are displayed for χ = 0.5, 2.5 and 10 in Fig 1C.For χ = 0.5, the plasmid slightly tends to be localized at the cell center, but it is still dominated by thermal fluctuations (Fig 1C1).The plasmid was stably localized at the cell center for χ = 2.5 (Fig 1C2), whereas for χ = 10, it showed directed movement and then reflected at the end walls, resulting in pole-to-pole oscillation (Fig 1C3).
Fig shows simulation results of Eqs 3 and 4 for χ = 2.5 (S3A Fig) and χ = 10 (S3B Fig) with K d = 0.001.In the former case, the plasmid maintained its location, whereas in the latter it traveled on the u(x) wave and moved unidirectionally.Fig 2A shows a steady velocity (= v) profile of plasmid movement for 0 < χ(= kN/(2l b )) < 10 and 0 < εð¼ DN=ð2l b ÞÞ < 10.In the case without thermal fluctuation, as χ and ε increased, the velocity monotonously increased above a critical curve on the χ − ε plane (Fig 2A).Above the curve, the plasmid showed directed motion.To analytically examine the change in the steady solutions of plasmid dynamics against χ values, we simplified Eqs 3 and 4 assuming u(x) � K d over x 2 [0, L] resulting in uðxÞ m ðK m d þuðxÞ m Þ ! 1 (S2 Text).Furthermore, by introducing a co-moving frame with a space-time coordinate z ≔ x − vt, where v is the steady velocity of the plasmid, defining u(x, t) ≔ U(x − vt, t), and solving the steady-state equation S3 and S4 Eqs in S2 Text, we obtained a relationship between v and χ (Fig 2B, green solid line) and the steady-state solution U st (z) (Fig 2C, green solid line) (S2 Text).Solutions for directed movement (|v| > 0) emerged at χ = χ c � 3.1 as a result of a (supercritical) pitchfork bifurcation, whereas the localized solution without motion (v = 0) existed over 0 � χ � 10.Hence, there were three solutions for χ > χ c (S2 Text and S4 Fig).Fig 2C shows steady solutions U st (z) for localization (purple) at χ = 2.5 and directed movement (green) at χ = 10.Note that the plasmid location is fixed at the origin (z = 0) on the space-time coordinates (S2 Text).For the case of v = 0, the shape of U st (z) was symmetrical and had its minimum at the origin, reaching an equilibrium state of the plasmid location (Fig 2C, purple).In contrast, the symmetry of U st (z) was broken for the case of |v| > 0, supporting a non-equilibrium traveling wave (Fig 2C, green).Interestingly, the minimum of the latter traveling wave was positioned at a location shifted from the origin of the plasmid (Fig 2C, inset figure, green), suggesting that the plasmid was "surfing" on the traveling wave.Next, we numerically calculated the steady velocity using Eqs 3 and 4 with K d = 0.001 over 0 � χ � 10 and confirmed the emergence at χ = χ c and the stability over χ > χ c of the solution of the equation S1 and S2 Eqs in S2 Text for the plasmid surfing on the traveling wave (Fig 2B, red dots).Therefore, these results demonstrated the existence of solutions for plasmid surfing on the traveling wave of ParA-ATP for a 1D case.

Fig 2 .
Fig 2. Chemophoresis engine mathematically validates plasmid surfing on the traveling wave of ParA-ATP.(A) Steady velocity (|v|) profile of plasmid movement for 0 < χ < 10 and 0 < � < 10 without thermal fluctuations.The plasmid starts moving above a critical curve on the χ − � plane.The white dotted line shows a parameter region in Fig 2B.The magenta dots at (χ, ε) = (2.5, 5), (10, 5) corresponds to parameter values for the steady solutions shown in Fig 2C.K d = 0.001, N = 40, and L = 40.(B) Relationship between v and χ in analytical (green for v 6 ¼ 0 and purple solid line for v = 0) for S3 and S4 Eqs in S2 Text, and simulated (red dots) solutions for Eqs 3 and 4. Solutions for directed movement (|v| > 0) emerge at χ = χ c � 3.1 as a result of a supercritical pitchfork bifurcation, whereas a solution for localization (v = 0) exists over 0 � χ � 10.K d = 0.001, N = 40, and L = 40.(C) Analytical solutions of the ParA-ATP pattern U st (z) in localization (purple) at χ = 2.5 and directed movement (green) at χ = 10 for S6 and S7 Eqs in S2 Text.The inset figure shows an enlarged view of U st (z) for z/L 2 [−0.02: 0.02].The plasmid location is fixed at the origin (z = 0) on the space-time coordinates.K d = 0.001, N = 40, and L = 40.(D) l b dependency of the directed motion of the plasmid for the analytical solution S7-S9 Eqs in S2 Text (solid lines) and numerical result calculated from Eqs 3 and 4 (black dots).For N = 20, 30, 40, 50, there exists a solution with the directed movement, whose velocity monotonously decreased with l b , and an inverse pitchfork bifurcation occurred at a critical value of l b , resulting in the only solution with v = 0.The numerical results showed slight deviation from the analytical solution for the ranges of small and large l b , suggesting the breakdown of the approximation u(x)�K d .For clarity, the numerical result for v = 0 was displayed only in the case of N = 40.k ¼ 0:1; D ¼ 0:05, and L = 40.https://doi.org/10.1371/journal.pcbi.1010324.g002 p.This analytical representation has a bifurcation at k ¼ k c ≔2=½NðDN=2 þ 1Þ�: for k larger than this critical value, there exist two realvalued solutions representing the directed movement.Corresponding to Fig 2A and 2B, we obtained the phase diagram and the relationship between |v| and k for D ¼ 0:05 and N = 40, for the parameter values adopted in Figs 1 and 2 (S6 Fig).The bifurcation point k c (= 0.025) in the vanishing size limit (S6B Fig) is smaller than the corresponding value (k c � 0.031) for l b = 0.2 in Fig 2B, suggesting that the directed motion for a larger l b requires a larger k, and probably N as well.We also examined l b dependency of the directed motion of the plasmid for the cases N = 10, 20, 30, 40, 50, using the analytical solution S7-S9 Eqs in S2 Text (Fig 2D).There were no solutions showing the directed movement for N = 10, whereas, for N = 20, 30, 40, 50, there exists a solution with the directed movement, whose velocity monotonously decreased with l b , and then an inverse pitchfork bifurcation occurred at a critical value of l b , resulting in the only solution with v = 0 (Fig 2D).Next, we numerically calculated l b dependency of the velocity and confirmed the disappearance of the traveling wave solution with the increase in l b (Fig 2D, 4ðmmÞ, while ρ is a dimensionless relative density normalized by the cytoplasmic one of ParA-ATP assigned to u 0 (See also S1 Text for normalized parameters and S1 Table for in vivo parameter values).We have performed simulations for the single plasmid (M = 1, S7 Fig upper) and the twoplasmid (M = 2, S7 Fig bottom) cases using Neumann boundary condition as in Fig 1.We investigated how the directed (or pole-to-pole oscillatory) movement of plasmids emerges with the increase in ρ, by calculating time-averaged values of velocity for different cell lengths (S7A Fig).The directed motion appeared at a ρ = ρ c for cell lengths greater than 0.8(μm) for M = 1 (S7A1 Fig) and 1.6(μm) for M = 2 (S7A2 Fig), Fig 2B and 2C) and the stability (S5 Fig) of the plasmid-surfing pattern.The solution emerged through the symmetry-breaking transition of the ParA-ATP spatiotemporal pattern at a critical χ.We mathematically showed the directed movement emerges even in the limiting case of vanished plasmid-size l b !0 (S6 Fig).Also, with an increase of the plasmid size l b , the solution for the directed movement disappeared as a result of an inverse pitchfork bifurcation (Fig 2D).By using parameters capturing in vivo conditions, we demonstrated that the chemophoresis engine can work even in bacterial cells (S7 Fig).The plasmid surfing also worked for a two-dimensional (2D) case (S8 and S9 Figs).The simulation results for the 2D case in Eqs 3 and 4 are shown for χ = 10 (S8A Fig) and χ = 50 (S8B Fig).In the former case, the cargo maintained its location, and u(r) had a symmetrical shape (S9A Fig), whereas directed motion by surfing on an asymmetrical traveling wave of u (r) was observed for the latter (S9B Fig), just like the 1D case.
and oscillatory waves (left).The oscillatory behavior of plasmids does not disrupt time-averaged equi-positioning.Steady multi-modal distributions of plasmids are sustained (Compare (B) right and (C) right).K d = 0.1, ε = 5, and L = 5.The distributions (right) were generated using 10 7 samples over 10 5 time step.(PDF) S3 Fig. Simulation results of Eqs 3 and 4. χ = 2.5 (A) and χ = 10 (B) with K d = 0.001.The red lines show the spatial pattern of u(x, t) (=ParA-ATP) for each t. Green dots show a plasmid location for each t.In the former case (A), the plasmid maintains its location, whereas it surfs on the traveling wave of u(x, t) and moves unidirectionally in the latter (B).These results were compared with the analytical results in the main text and in Fig 2. ε = 5 and L = 40.(PDF) S4 Fig. Analytical solution of εΔμ(v, χ) − v for visualization of the self-consistent equation S9 Eq in S2 Text.S9 Eq in S2 Text shows a pitch-fork bifurcation at χ = χ c � 3.1, and has three solutions for χ > χ c .ε = 5 and L = 40.(PDF) S5 Fig. Linear stability analysis of the surfing-on-wave pattern U st (z) and z st ξ ð¼0Þ against external disturbances.