Engineering transient dynamics of artificial cells by stochastic distribution of enzymes

Random fluctuations are inherent to all complex molecular systems. Although nature has evolved mechanisms to control stochastic events to achieve the desired biological output, reproducing this in synthetic systems represents a significant challenge. Here we present an artificial platform that enables us to exploit stochasticity to direct motile behavior. We found that enzymes, when confined to the fluidic polymer membrane of a core-shell coacervate, were distributed stochastically in time and space. This resulted in a transient, asymmetric configuration of propulsive units, which imparted motility to such coacervates in presence of substrate. This mechanism was confirmed by stochastic modelling and simulations in silico. Furthermore, we showed that a deeper understanding of the mechanism of stochasticity could be utilized to modulate the motion output. Conceptually, this work represents a leap in design philosophy in the construction of synthetic systems with life-like behaviors.

chains. 1 g of Phe-terminated copolymer was weighed into a Schlenk flask and dissolved with ca. 3 mL of dry DMF and cooled in an ice bath. To the cooled mixture 160 mg of NCA-BLG was added under Ar and the reaction was left under a constant flow of N2 for 24 h. The product was precipitated into cold methanol and analysed by 1 H NMR to confirm the overall composition and, in particular, the presence of benzylic and aromatic protons at 5.0-5.2 and 7.1-7.4 ppm, respectively. Benzyl-protected terpolymer was dissolved in 10 mL THF and 10 mL methanol was added before applying hydrogenation using the H-Cube at 60 °C with 30 bar of H2 pressure and a flow rate of 1 mL min -1 to facilitate removal of the benzyl protecting groups. The product was concentrated, then precipitated into cold ether, and dissolved in dioxane before lyophilisation to yield a waxy solid, 0.9 g (85 % yield). 1 H NMR was used to confirm successful deprotection of the PBLG units, and GPC data indicated that the polydispersity did not increase beyond 1.1 during this process.
Synthesis of azido-poly(ethylene glycol)-b-poly(ε-caprolactone)-g-poly(trimethylene carbonate) (azido-PEG68-b-PCL50-g-PTMC50) block copolymers (azido-pol) The synthesis protocol was modified from a literature procedure.2 First, azido-PEG-OH macroinitiator (3 kDa, 0.1 mmol, 300 mg) was weighed into a round bottom flask along ε-CL (5 mmol, 553 μL) and TMC (5 mmol, 510 mg) and dried via azeotropic evaporation of added toluene (x3). The dried reagents were then re-dissolved in dry DCM (10 mL) and dry methanesulfonic acid (0.2 mmol, 14 μL) was added under argon. The reaction was stirred at 30 °C overnight, after which there was no sign of unreacted monomers confirmed by 1H NMR. After the reaction was finished, the solution was concentrated and copolymer precipitated into ice cold MeOH (x2), followed by dissolution in dioxane and lyophilization to yield 1.01 g of a waxy solid (78% yield). Copolymer composition was ascertained from the 1H NMR spectrum in the same manner as for step 1 of the terpolymer synthesis above. GPC analysis (using a PL gel 5 µm mixed D column, with THF as eluent and PS standards) yielded Mn = 15.7 kDa, Đ = 1.12. Copolymer composition was confirmed by 1H NMR (Supplementary Fig. 4).

Enzyme modification and activity test
Catalase modification: synthesis of mCAT(Cy5-DBCO-Catalase) Catalase was modified with Cy5 and DBCO as follows: 20 mg catalase was dissolved in 2.5 mL 0.1 M sodium bicarbonate buffer (adjusted to pH = 8.2) in a glass vial, to which 125 µL of an aqueous solution of sulfo-Cy5-NHS ester (4 mg/mL) was added. This 8:1 (dye: protein) stoichiometry was chosen to ensure sufficient dye labeling as the NHS ester is easily hydrolyzed. The mixture was covered with aluminum foil and stirred overnight at 4°C and then dialyzed for 24 hours against 1x PBS buffer. Cy5 labeled catalase was purified via fast protein liquid chromatography (FPLC, BioRad NGC system). DBCO-PEG4-mal was dissolved in PBS by vigorous vortexing and sonication and added into Cy5-catalase solution. The reaction was carried out overnight by gentle stirring at 4°C. The obtained mCAT was purified by FPLC. The concentration was determined to be 2.3 mg/mL by Nanodrop. Cy5 labelling number per catalase tetramer was determined to be 2 by UV-vis spectroscopy. The presence of DBCO on catalase was confirmed by fluorescence emission after addition of the fluorogenic agent 3-azido-7hydroxycoumarin. 3-azido-7-hydroxycoumarin is highly quenched and upon azido-alkyne reaction with DBCO, the fluorescence is restored.
Urease modification: synthesis of mUR(Cy5-DBCO-Urease) Urease was modified with Sulfo-Cy5-NHS ester and DBCO-NHS ester as follows: 100 mg urease was dissolved in 10 mL 0.1 M sodium bicarbonate buffer (adjusted to pH = 8.2) in a glass vial, to which 1.14 mg sulfo-Cy5-NHS ester powder and 0.3 mg DBCO-NHS ester (in DMSO)were added. The excess of NHS ester was to ensure sufficient dye labeling as NHS ester is easily hydrolyzed. The mixture was covered with aluminum foil and stirred overnight at 4°C and then dialyzed for 24 hours against 1x PBS buffer. The obtained mUR was then purified with FPLC. The concentration was determined to be 3.7 mg/mL by Nanodrop. Cy5 labelling number per urease was determined to be 1.4 by UV-Vis spectroscopy. The presence of DBCO on mUR was confirmed by fluorescence emission after addition of the fluorogenic agent 3-azido-7hydroxycoumarin.
Bovine serum albumin modification Bovine serum albumin (BSA) was succinylated and modified with dye Atto 488 as follows: 5 mg BSA was dissolved in 2.5 mL 100 mM NaHCO3 (adjusted to pH 9) in a 1.5 mL Eppendorf tube (BSA concentration 2.5 mg/mL), to which 20 μL Atto 488 -NHS ester (10 mM, in DMSO) was subsequently added. The reaction mixture was gently stirred at room temperature overnight. The next day, succinic anhydride was added directly to the reaction mixture at a final concentration of 50 mM. The reaction mixture was gently stirred under aluminum foil for another 4 hours before purification by FPLC and concentration to approximately 1.2 mg/mL. Atto 488 labeling number per BSA was determined to be 1.2 by UV-Vis spectroscopy.
Catalase activity assay before and after modification The activity of mCAT was determined by measuring the decomposition of hydrogen peroxide over time via UV-Vis spectroscopy (JASCO V-650). In these assays, 0.1 mL of diluted catalase (unmodified or modified) in PBS solution was mixed with 2.9 mL hydrogen peroxide solution (0.036 wt% in water), and the characteristic absorbance at 240 nm (A240) was monitored over time. Enzyme activity was calculated according to: where 3.45 corresponds to the decomposition of 3.45 µmoles of hydrogen peroxide in a 3.0 ml reaction mixture producing a decrease in the A240 from 0.45 to 0.40, df is the dilution factor and time (minutes is the time for A240 to decrease from 0.45 to 0.40. Enzyme activities of unmodified catalase and modified catalase were determined after 3 individual measurements to be about 10×10 3 U/mg and 5×10 3 U/mg respectively.
Urease activity assay before and after modification Urease activity was tested by measuring pH change in solution upon addition of its substrate urea. First the pH of urease solution was adjusted to around 4.0 by HCl addition, then urea was added and pH change over time was recorded. Ten minutes after urea addition, the pH recording stopped. Urease activity before and after modification were both tested. After modification, urease kept about 74% activity according to: Relative activity = ∆pH (urease after modification) ∆pH (urease before modification) .
Enzyme activity test after surface attachment via SPAAC reaction Catalase activity was tested after surface attachment qualitatively. After surface attachment and removal of unbound catalase via centrifugation, 0.2 mL 5 wt% H2O2 was added to 0.2 mL mCAT-coacervate suspension in PBS buffer, and the generation of a large amount of oxygen bubbles confirmed catalase was still active after surface attachment.

Fabrication of mCAT or mUR functionalized coacervates
Coacervates were prepared based on a sonication method as follows: to 200 µL 0.5 mg/mL quaternized amylose (Q-Am) solution in an Eppendorf tube was added 100 µL 0.5 mg/mL carboxymethyl amylose (Cm-Am) solution to induce coacervation. The mixture was placed in a sonication bath for 2 minutes to allow the coacervate growth, followed by an addition of 10 µL 50 mg/mL polymer mixture (20%v/v 50 mg/mL azido-PEG68-b-PCL50-g-PTMC50 and 80%v/v 50 mg/mL terpolymer in PEG 350) to generate a polymeric membrane and stabilize the coacervate core. The size of such coacervates was effectively limited to 1-3 µm by utilizing sonication during coacervate droplet growth. Coacervates larger than 3 µm could also be generated by using a method based on employing an Eppendorf shaker as previously reported. 1 Large coacervates were prepared via the sonication method stated above, except the amylose mixture was placed in an Eppendorf shaker instead of a sonication bath. A broad range of coacervate sizes was achieved through the combination of changing mixing methods (shaker or sonication), altering volume and growth time of the coacervate core. Supplementary Figure 14 shows the size distribution of coacervates with diameters of 1.2 ± 0.4 µm, 2.6 ± 0.7 µm, 3.6 ± 1.6 µm and 10.3 ± 3.4 µm.
After the assembly of coacervates, certain amounts (listed below) of mCAT (Cy5-DBCO-Catalase) or mUR (Cy5-DBCO-Urease) were added to allow mCAT / mUR attachment on the polymeric membrane through strain-promoted azide-alkyne cycloaddition (SPAAC) between the azido-group installed on the polymer membrane and the DBCO moiety introduced in the enzyme.
The reaction was carried out for two hours, followed by separation of unbound enzymes by centrifugation and refreshing the supernatant with buffer. Surface enzyme density (Fig. 4) was tuned by changing the added amount of mCAT or mUR. For mCAT-coacervates, low, medium and high surface density corresponded to 4.6, 46, or 81 µg mCAT addition. For mUR-coacervates, low, medium and high surface density corresponded to 5.4, 54, or 108 µg mUR addition. The actual density of both enzymes on the surface of coacervates was calculated using fluorescence spectroscopy (vide infra).

Size analysis of coacervates
The size of the coacervates was obtained by confocal image analysis. To facilitate size analysis, succinylated Atto 488 modified BSA (sBSA) was loaded in the interior of the coacervates. To assemble sBSA loaded coacervates, after addition of Cm-Am, 10 µL sBSA was added to the amylose mixture, followed by a growth time of 2 min and addition of the membraneforming terpolymer.
Confocal laser scanning microscopy (Leica TCS SP8) was used to capture images of sBSAloaded coacervates with a 488nm laser line using a 63×, 1.20 NA water immersion objective. The pinhole was set to 1 Airy Unit (156 μm). The images were then converted to binary, and the area of each coacervate was analyzed using software ImageJ. The diameter of coacervates was obtained by averaging at least 100 coacervates.

Motility test for mCAT-/mUR-coacervates
Experimental chamber A simple experimental chamber was designed and prepared to minimize side effects that could be mistaken as self-propulsion, such as drift or solution evaporation. This chamber was made from two glass microscopy slides spaced by two pieces of autoclave tape. Autoclave tape was first attached to a larger glass slide on two ends, followed by the addition of sample in the middle, and capping by a smaller glass slide.

Optical recording
The videos of coacervate motion were recorded using a camera (Hamamatsu Digital Camera C11440) on an inverted optical microscope (Leica DMi8). A 63× water immersion objective was used for this recording. The coacervate suspension and their corresponding substrate (hydrogen peroxide or urea) were first mixed thoroughly and immediately added to the chamber. The recording started shortly after. For mCAT-coacervates, hydrogen peroxide concentration after mixing was set to 10 mM. This very low concentration was chosen to avoid oxygen bubble generation that could disturb the liquid inside the chamber and mask self-propulsion. Sample in the experimental chamber was replaced at least every 3 min to maintain the same substrate concentration.

Data analysis of motion
A tailor-made Python script was used to track the coacervates and obtain the mean square displacement (MSD). MSD is a measure of deviation of the position of a particle with respect to its initial position over time. It is commonly used to analyze the dynamics of self-propelled particles and is calculated as below in 2D projection: MSD ( ) = 〈( ⃗( ) − ⃗ (0)) 〉, where ⃗(0) is the initial position of the coacervate, and ⃗( ) is the position of the coacervate when time is t. The obtained MSD was plotted against the time interval Δt. To extract the velocity of the particles in the self-propulsive regime, a least-squares fitting to the equation MSD(Δ ) = 4 ∆ + ∆ was performed using OriginLab software.
In situ crosslinking of mCAT on the polymeric membrane After mCAT coacervate fabrication, 109 µL of glutaraldehyde dissolved in 1x PBS (GA, 275 mM) and 91 µL of 1x PBS were added to 100 µL mCAT coacervate suspension. The final concentration of GA was 100 mM. The reaction was carried out for 3 hours at room temperature with mild stirring. Unreacted GA was removed by centrifugation and resuspension twice.
Fluorescence recovery after photobleaching of enzyme-attached coacervate membrane To monitor the lateral diffusion of surface-attached enzymes (mCAT and mUR), fluorescence recovery after photobleaching (FRAP) (Leica TCS SP8 confocal microscope) was performed on the membrane. A selected circular area on the membrane was bleached (100% laser intensity), and the fluorescence intensity in the bleached area was monitored (1-2% laser intensity) over time. Faster fluorescence recovery in the bleached area indicated faster lateral diffusion of enzymes on the membrane. As micron-sized particles have a characteristic rotational time (the time required for a particle to rotate 360° in solution around its axis) due to Brownian motion, we chose coacervates with diameters of 10 µm (with a characteristic rotational time of ~500 s) for the FRAP measurements (substantial recovery observed after 52 s in Fig. 2D) to minimize the effect of rotation on apparent lateral diffusivity.
Data was analyzed as follows: Fluorescence intensity at time t, I(t) was background subtracted and corrected for unintentional photobleaching, which was calculated as where ROI(t) is the average intensity of the bleached area at time t, Ref(t) is the average intensity of an unbleached fluorescent area (same size as bleached area) and Bg is the average intensity of a background area. I(t) was further normalized such that pre-bleach intensity was set to 1. This was done by dividing I(t) by average background-subtracted pre-bleach intensity within the bleach area.
Normalized intensity I(t)normalized was plotted against time to get the recovery curve. This recovery curve was fitted with I(t)normalized = B -Ae -τt , assuming exponential kinetics, to obtain parameter τ. In the above equation, τ is the recovery time constant and A and B are two constants. Then recovery half-time τ1/2 and apparent diffusion coefficient DL were calculated from τ1/2 = -ln0.5 / τ and DL = 0.88 ω 2 / 2τ1/2 (for a circular bleaching area) respectively. 4 We therefore obtained a DL of 0.0352 ± 0.0008 µm 2 / s for mCAT-coacervates and a DL of 0.0298 ± 0.0005 µm 2 / s for mUR-coacervates. The results were obtained by measuring at least 6 different enzyme-attached coacervates for both cases.
FRAP was performed on crosslinked mCAT-coacervates following the same protocol. The recovery curve was plotted in the same way, however, the fitting equation for non-crosslinked mCAT coacervate was not suitable for the crosslinked coacervates because the interaction between catalase enzymes (because of crosslinking) complicated the fluorescence recovery, which was not taken into account in the fitting equation  Synthesis of PEG44-b-PCL50-g-PTMC50-b-PGA8 terpolymer. Poly(ethylene glycol)44 monomethyl ether was used to initiate the ring opening polymerisation of ε-caprolactone and trimethylene carbonate (step 1). The terminal alcohol of this polymer was subsequently modified via a Steglich esterification with Boc-L-Phe-OH to yield a primary amine after TFA deprotection (step 2). The final poly(L-glutamic acid) block was introduced by the ring opening polymerization of Ncarboxyanhydride γ-benzyl-L-glutamate, followed by hydrogenation (step 3).   The formation of terpolymer stabilized coacervates proceeds via hierarchical self-assembly. First, oppositely-charged modified amylose polymers are mixed, immediately forming complex coacervate microdroplets that coalesce and increase in size over time. This growth is arrested by the addition of a carefully-designed block terpolymer, PEG-PCL-PTMC-PGlu, which selfassembles on the surface of coacervate microdroplets (A), driven by electrostatic attraction between PGlu (negative) and Q-Am (positive). Large, stable populations (B) of coacervates are thus obtained.   FRAP measurements were performed on mUR-coacervates to test the lateral diffusivity of urease. (A) Substantial recovery 52.2 seconds after laser bleaching was observed (Green: Cy5 from mCAT). (B) FRAP recovery curve with black dots representing experimental data (n = 5, mean ± SEM) and red curve exponential fitting. Lateral diffusivity of surface-bound mUR was determined to be around 0.03 µm 2 /s. Scale bar represents 2 µm. Three individual batches of coacervates were prepared for FRAP measurements.   Schematic illustration of relationship between vector density (enzyme density in 2D) and total generated propulsion. Please note: this figure is only for illustrating the idea of interplay, and it doesn't imply the moving direction of enzymes or coacervates.   The size of coacervates can be finely tuned by adjusting amylose derivative volume and altering the preparation method (shaker or sonication). Coacervates with diameter of 1.17 ± 0.39 µm were prepared by the sonication method with growth time of 1.5 min; coacervates of 2.58 ± 0.71 µm were prepared by the shaker method with 0.2 unit volume (40 µL Q-Am + 20 µL Cm-Am); coacervates of 3.59 ± 1.58 µm were prepared by the shaker method with 0.4 unit volume (100 µL Q-Am + 50 µL Cm-Am); coacervates of 10.3 ± 3.40 µm were prepared by the shaker method with 1 unit volume (200 µL Q-Am + 100 µL Cm-Am). For all cases, n > 100. Three individual batches of coacervates were prepared for size measurements. Scale bar represents 5 µm.    We consider a spherical particle of radius . Enzymes are attached on the surface of the sphere and can diffuse along the surface with lateral diffusion coefficient . Due to the Brownian motion of the enzymes at the surface, the distribution of enzymes fluctuates around the homogeneous equilibrium concentration .
We assume that the particle velocity, is proportional to the instantaneous dipolar distribution of the enzymes: (1) The proportionality constant is a characteristic velocity scale per enzyme that depends on the specific phoretic propulsion mechanism, represents the total number of enzymes and the unit vectors , with = , , , define the reference frame of the particle and , with = , , , define the components of the dipole. The dimensionless vector ( ) defines the instantaneous direction and magnitude of the dipolar distribution relative to the equilibrium concentration Γ . The dynamics of the particle position vector are given by the superposition of the selfpropulsive velocity and the translational Brownian motion: With the translational noise due to Brownian motion given by a Gaussian random variable with zero mean and delta correlated in time 〈 ( ) (´)〉 = 2 ( −´), with the translational diffusion coefficient of the particle. This model is reminiscent of the Active Brownian Particle (ABP) model, with the difference that the dipole ( ) changes due to the surface diffusion of the enzymes and due to the rigid-body rotations of the particle.
To close the problem, we need to evaluate the evolution of the dipolar concentration of the enzyme ( ). To do so we derive a Langevin equation for the vector ( ).

Langevin equation for the fluctuating dipole
We assume that the fluctuations of the enzyme concentration are small and can be computed as if the system was in equilibrium. 9 Considering the full non-equilibrium features of the fluctuations would couple different spatiotemporal modes and different fields, which would make the problem extremely complicated. We believe that these interesting aspects deserve a separate analysis.
Using the framework of linear fluctuating hydrodynamics, it can be shown that the dipolar concentration vector, scaled by equilibrium distribution , obeys the Langevin equation where the noise ( ) represents the change of dipole orientation due to the rigid-body rotational Brownian motion of the particle. This noise is Gaussian, it has zero average and it is deltacorrelated in time 10,11 where is the rotational diffusion coefficient of the particle, which is given by the Stokes-Einstein formula where represents the thermal energy and is the viscosity of the liquid suspending the particle The surface diffusion of the enzymes contributes to the fluctuations of the dipole through the term ( ), which is Gaussian-distributed with zero mean and delta-correlated in time The first term on the right-hand side describes the tendency of entropy to restore a homogeneous distribution of enzymes. We The symmetric tensor − ( ) ( ) selects the components of the noise that are perpendicular to With the noise having variance given by The above equation states that the dipole orientation fluctuations are due to the sum of the rigid body rotation of the sphere and the surface diffusion of the enzymes. The characteristic time associated to the change of orientation of the instantaneous dipole is given by The correlation of the orientation at different times is therefore given by 10 where we denoted the modulus of ( ) as ( ) and we used the fact that The noise is a scalar and its variance can be shown to be given by The above equation states that the modulus of the dipolar enzymatic concentration is given by an Ornstein-Uhlenbeck process. 11 Straightforward integration of the equation yields the correlation 11 with the characteristic time = . The same-time magnitude variance is inversely proportional to the total number of enzymes at the surface . Finally, it can be shown that the modulus and the orientation of the dipole are two uncorrelated random processes.

Evaluation of MSD
We rewrite the equation of motion of the particle as ( ) = ( ) ( ) + ( ).
The solution of the equation is given by With ( ) is an increment of a Wiener process with zero mean and variance 〈 ( ) (´)〉 = 2 ( −´). The mean square displacement (MSD) of the particle is calculated as By plugging in the definition of the correlations of ( ), ( ) and of the Wiener process ( ) we obtain The MSD has the same functional form as that found in the case of the Active Brownian Particle (ABP) model. 12 Note that this relation is equivalent to Equation (1) in the main text and Equation (14) in Supplementary Note 2. We find a ballistic regime at small times and a diffusive regime at long times. The main difference between the MSD predicted by the ABP and that predicted by the equation above is given by the characteristic timescale over which the particle transitions from the ballistic regime to the diffusive regime, * = , is set by the inverse an effective rotational diffusion coefficient, given by While for the ABP model this characteristic timescale is simple given by . Compared to the ABP model, * depends on the surface diffusion coefficient and on the total number of adsorbed enzymes . Usually, there are many enzymes adsorbed on the surface of the particle and ≪ 1.
In the case of mobile enzymes, their surface diffusion contributes to the change of the dipole orientation, which determines the direction of the self-propulsion. This means that there are two mechanisms responsible for the change of the dipolar distribution of enzymes: (i) diffusion of the enzymes along the surface and (ii) rigid-body reorientation of the entire particle. These two mechanisms act in parallel and either one or the other can dominate, depending on the parameters.
In the case 1 + ≫ , the dipolar fluctuations of the enzyme distribution due to surface diffusion are faster than the rigid-body rotational diffusion and they determine the timescale over which the particle transitions from ballistic motion to diffusive motion * ≈ . Conversely, in the case ≫ 1 + , the enzymes' diffusion along the surface is slower than the rigid-body rotational diffusion of the particle. In this case, the reorientation of the dipole is mainly due to rigid-body rotations and the characteristic timescale recovers that predicted by the ABP model * ≈ .
Finally, since the experiments are performed in two dimensions, we restrict the MSD obtained above to the case of motion restricted on a plane. In this case, the MSD obtained above is multiplied by a factor 2/3, see Equation (7) of the next section, which yields: This equation can be readily used to fit the mean square displacement measured in the experiments. It can be further simplified for 2-parameter fitting procedure: where DT can be obtained from the Stokes-Einstein formula: We obtained A = 0.84 ± 0.14 and Deff -1 = 1.00 ± 0.12 for fitting shown in Supplementary  Figure 12.
In the next section, numerical simulations of the Brownian motion of each enzyme along the surface of the active particle are presented. The velocity of the particle is calculated from the diffusion of each individual enzyme, which allows to confirm the analytical theory presented here and to explore cases where the number of enzymes is small and fluctuations are large.

Supplementary Note 2: Stochastic model for motility of a protocell with dynamic membrane
In Appendix S1 we derived an analytical expression for the mean squared displacement (MSD) profile over time. The reader should realize that an individual's particle squared displacement profile will deviate from this mean as a result of stochastic variability. In this Appendix we present a stochastic simulation for the motility of an individual protocell with dynamic membrane. The analytical expression for the MSD can then be verified by comparison to the average of simulated squared-displacement profiles.

Spherical Brownian motion
Simulation of a spherical Brownian motion (SBM) is key in this work. A SBM starting in y ∈ S 2 (R), can be defined by the diffusion equation of the transition density x → ρ y (x, t) as where ∇ 2 S 2 (R) is the Laplace-Beltrami operator on the sphere S 2 (R) and D > 0 the diffusion coefficient. We have adopted the simulation strategy presented by Mijatović et al. (2020), 13 but instead we used a normal approximation to simulate from the Wright-Fisher distribution as discussed by Jenkins and Spano (2017). 14 The latter approximation is appropriate since we will simulate for small time steps. Simulation of the SBM for a time step ∆t is performed by Algorithm 1. The mean distance traveled from the starting points is 0 as illustrated in Supplementary Fig. 18. The mean squared distance from the starting point does depend on the diffusion coefficient and is illustrated in Supplementary Fig. 18 for D ∈ {0.03, 0.1, 0.5}. The signed distance from the starting point is computed as R · ∆σ · sign(φ 2 − φ 1 ), where ∆σ is the angle between the vector to the starting point and the vector to the new point, φ 1 and φ 2 equal the azimuthal angle of the spherical coordinates of the starting point and the new point respectively. The ∆σ is known as the central angle and can be computed as where λ 1 and λ 2 equal the polar angle of the spherical coordinates of the starting point and the new point respectively. Figure 18: The expected (signed) distance from the starting point (left) and the expected squared distance from the starting point, for diffusion on a sphere with R = 5 (middle) or R = 1 (right) and D ∈ {0.03, 0.1, 0.5}, based on 10000 simulations.

Motile particle ignoring rotational and lateral diffusion
Let us consider the two-dimensional stochastic process X t that represents the position of a motile particle at time t in the (x, y)-plane, assuming that there is no rotational diffusion of the particle, nor lateral diffusion of the enzymes attached, i.e., where B x (t) and B y (t) independent one-dimensional Brownian motions and D T represent the translational diffusion coefficient. 12 These Brownian motions model the translational diffusion of the particle. Furthermore, where N equals the random number of enzymes attached to the particle, x i equals the position vector of enzyme i projected on the z = 0 plane and v c is a normalization constant that could change with the presence of substrate (fuel). Note that v0 xi R then represents the effective-propulsion-direction vector. It is important to realize that it is unknown whether v c should be negative (motor pushes the particle ) or positive (motor pulls the particle). However, both scenarios will result in the same mean squared displacement (MSD) curve.
The MSD, in the (x, y)-plane, over time of such particles equals Since x i and x j , for i = j, are independent and E[x i ] = 0 due to the symmetry of the particle. By Wald's identity, the expectation of v 2 0 can be written as where κ = E |x i | 2 and λ = E [N ]. The projection of velocity vector of the enzyme attached in coordinate (x, y, z) on the sphere to the z = 0 plane equals (x, y). Under the assumption that the particle is functionalized randomly, we know that Thus, We have simulated 1000 times X t for 10 seconds, where we have assumed R = 1, v c = 1, √ 2D T = 1, N ∼ P oi(λ) for λ ∈ {20, 100, 1000}. The MSD in the (x, y)-plane of these simulations, as well as the theoretical curve in (30), are presented in Supplementary Fig. 19 and do match. It is important to remark that the MSD curve is not unique if D T is given and depends on the product of v 2 c and [N ], i.e. lowering the expected enzyme density with a factor c results in the same MSD curve if v c is increased by a factor √ c.
Figure 19: MSD for a particle ignoring lateral diffusion of the enzymes and rotational diffusion of the particle. Curves are presented for λ ∈ {20, 100, 1000} (black), where the MSD increases with the expected enzyme density. Theoretical curves are shown as dotted gray lines and nicely overlay with the means derived from the simulation (n = 1000).

D T estimates
The diffusion of small (R = 0.6um) and large (R = 1.8um) particles has been measured experimentally in the absence of fuel (see Materials). For those particles the MSD equals 4D T t, which allows us to estimate D T,0.6 and D T,1.8 . By fitting a simple (no intercept) linear regression we estimated translational diffusion coefficients equal D T,0.6 = 0.236 (0.0015) and D T,1.8 = 0.079 (0.0018) respectively. These diffusion coefficients are known to equal where k B T is the thermal energy and η is the viscosity of the medium in which the particle diffuses. 12 All experiments have been performed in similar conditions such that we can assume kBT 6πη to be constant. Based on the D T,0.6 and D T,1.8 values we estimate this constant as 0.236 · 0.6 = 0.141 and 0.079 · 1.8 = 0.142 respectively. Therefore we will assume that kBT 6πη = 0.14 in the rest of this study.

Motile particle with lateral diffusive enzymes
We further extend the model to describe the stochastic process Y t that represents the position of a motile particle at time t in the (x, y)-plane, where enzymes (motors) additionally undergo lateral diffusion, and is modeled as where dt equals the size of the time steps of the random walks and v j·dt is the net-propulsion vector, in the (x, y)-plane, as a result of the enzymes at time j · dt such that v x,j·dt = N i=1 x i,j·dt . This net propulsion is not constant as a result of the lateral movement of the enzymes. We assume each enzyme to undergo an independent SBM with diffusion coefficient D L over the surface of the sphere (see Section 1 for a description). Notice that the speed of the diffusion is independent of the size of the particle; two examples of the movement of an enzyme on a particle with R = 0.6 and R = 1.8 respectively, for D L = 0.03 are illustrated in Supplementary Fig. 20. Figure 20: Diffusion of an enzyme during 10 seconds and D L = 0.03 on a spherical particle with R = 0.6 (left) and on a particle with R = 1.8 (right).

FRAP experiment
To validate how accurate this SBM describes the lateral diffusion of the enzymes we have mimicked the fluorescence recovery after photobleaching (FRAP) experiment that was performed in the lab. Firstly, we have randomly 'functionalized' a spherical particle with R = 5.0 using a Poisson-distributed (λ = 10000) number of enzymes N . Subsequently, part of the enzymes are hit with a 'laser' with radius 1.25nm centered around ( √ 5 2 − 1.25 2 sin(α), √ 5 2 − 1.25 2 cos(α)), where α = acos From the FRAP experiment performed in the lab it was observed that roughly a fraction of 35% of the enzymes within the x − y area that the laser hit were not photo-bleached. Therefore, we have independently removed each enzyme within the (x − y) laser area with a probability 0.65, see Supplementary Fig. 21. Then we count the variable number (due to the lateral diffusion) of enzymes within the bleached region, and within 1 micron of the focal plane z ∈ (−0.5, 0.5), over time, and compared this number to a reference region of the same size at the other side of the particle. In Supplementary Fig. 21 we also show the situation after 5 seconds. Figure 21: Particle direct after photo-bleaching (t = 0, left, where we observe a bright 'hole' in the top left side) and particle after five seconds of enzyme diffusion (t = 5, right, where the 'hole' has diffused away) From the experiment D L was estimated as 0.035 for the catalase (mCAT) enzymes. The intensity curve obtained from simulation with this particular value of D L was compared to the observed intensity curves from the real FRAP experiments as shown in Supplementary Fig. 22. We conclude that the SBM model can be used accurately to describe the lateral diffusion of enzymes. The mean squared displacement from the starting point of an enzyme for t ∈ [0, 10] with D L = 0.035 was presented in Supplementary Fig. 18.

MSD profile
Although the position vectors of one specific enzyme are strongly correlated over time, the position vectors of different enzymes are still independent over time. As a result the MSD can be written as In Supplementary Fig. 23 we present the MSD profiles from simulations of particles with R ∈ {0.6, 1, 1.8}, λ = 100 and D L ∈ {0, 0.03, 0.1, 0.25, 0.9, 3}. In this work (see Appendix S1) we derived an analytical form for (34), where D eff = D R + (1 + 2 λ ) D L R 2 and D R equals the rotational diffusion coefficient which will be discussed in Section 4. For D R = 0, Supplementary Fig. 23 confirms the validity of this expression in our simulation study. For the experimental particles studied λ > 5000 as a result of which D eff can be approximated by D R + D L R 2 . Then, the MSD profile does only depend on λ via v c √ λ, for which reason we can simulate with a lower number of λ and a higher value of v c to seriously reduce the computation time.

Motile particles with lateral diffusive enzymes and rotational diffusion
Finally, we extend our model for the motile particle by accounting for the rotational diffusion. The rotational diffusion is modeled by rotation of the coordinate system. The rotation of the coordinate system follows the rotation of the vector (0, 0, 1) that undergoes a SBM on the unit sphere characterized with D = D R , where D R = kBT 8πη 1 R 3 , referred to as b t . Subsequently, the net-velocity vector at time t, is rotated following this SBM as described by Algorithm 2.
Algorithm 2 Simulating rotational diffusion with diffusion coefficient D. 1: Simulate b t according to Algorithm 1 using D = D R . 2: Let v t be equal to the net-velocity vector of the particle without rotational diffusion and let |v t | = n t . 3: Let z t be equal to the unit directional vector such that n t z t = v t . 4: Set O(z t ) = I − 2uu , where u = ((0,0,1) −zt) |((0,0,1) −z)| . 5: As D R = kBT 8πη 1 R 3 , large particles rotate slower than small particles. 12 Due to the rotations the positions of the different individual enzymes over time become correlated. The MSD can be written as and the analytical form has been presented as Equations(35). In Supplementary Fig. 24 we present the influence of D R on the the MSD presented in (36), for particle with R ∈ {0.6, 1, 1.8}, where λ = 100 (and √ 2D T = 1), both for D L = 0.03 and D L = 0. Again we have validated the accuracy of the analytical MSD expression (35), presented as the green lines in Supplementary Fig. 24.
The theoretical expression (37) does indeed fit the MSD curves derived with our SBM model for the rotation diffusion, from which we conclude that the latter was indeed an appropriate approach.

Experimental data fitting and simulation
Based on the no-fuel data we did derive kBT 6πη as 0.14 as a result of which D T,0.6 = 0.233 and D T,1.8 = 0.078. The rotational diffusion coefficient is known to equal kBT where X = {low-mCAT, medium-mCAT, high-mCAT, low-mUR, medium-mUR, high-mUR} and Y = {mCAT, mUR}. The fitting is performed based on non-linear least-squares estimation using the nlsLM() function, from the package minpack.lm in [R]. Based on the FRAP experiments we use the restrictions D L,0.6,mCAT ≥ 0.036 and D L,0.6,mUR ≥ 0.03. The parameter estimates, standard errors and predicted MSD curves can be found in Supplementary Fig. 26

Simulation
Using the parameters obtained from the experimental data we can now simulate all types of particles that have been studied. The MSD curves of the simulations for the different settings are presented together with the experimentally observed MSD curves in Figures 27, 28 and 29. It is important to remark that with 10000 simulations the mean of the simulated particles is equal to the analytical fit presented in Supplementary Fig.  26.   Thus, it becomes clear that we can appropriately describe the experimental observations with the model proposed. To do so, the lateral diffusion coefficient in the case of the large mCAT particles should substantially differ from the lateral diffusion coefficient in the case of the small particles.

Instantaneous velocity
For the simulated particles from Section 5.1 the 'instantaneous' velocity was derived over time as The distribution of these instantaneous velocities are presented as Supplementary Fig. 30. The mean velocity of the simulated large particles (1.8 um/s) deviates from the mean velocity of the large experimental (non-crosslinked) particles (3.24 um/s). The velocity highly depends on the size of the translational diffusion coefficient D T . On the contrary, the influence of D T on the MSD was only minor for the large particles, such that the fits to the data (with fuel, Supplementary Fig. 27) were insensitive to the value of D T . In Supplementary Fig. 31 the cumulative probability distribution of the experimentally observed velocities is presented together with the distribution of simulated large particles using kBT 6πη = 0.45, V 1.8 = 0.576 and D L,1.8 = 0.917. For this simulation setting the mean velocity (3.18 um/s) corresponds to the experimentally observed average. It is important to note that the experimentally observed velocities are discretized as a result of the finite pixel size used during the imaging. We conclude that the variability in the velocities is well explained by the stochastic model. However, it becomes clear that a fraction of approximately 1% of the observed velocities is too extreme (> 12.5 um/s) to be explained by our stochastic model, as these velocities are never observed in the simulations. Figure 31: Cumulative distribution function, i.e. fraction of the time that the velocity is less than or equal to x, of the instantaneous velocity (left) observed experimentally (red) and derived from 10000 simulated (10s) large particles (black). Adjusted x-and y-axis are presented to emphasize the agreement (middle) and difference (right) of the two distributions.
In summary, we can conclude that we can explain the motile behavior of the vast majority of the particles. In small-sized particles the velocities are well described by the stochastic model. For large particles, 10% of the large coacervates modified with catalase show a run and tumble behaviour (Movie S6) that far exceeds the velocity distribution and cannot be explained by this model. We assume that a more complex combination of stochastic processes is at the core of this yet unexplained variability; for example due to interactions between the enzymes.