Passive and active colloidal chemotaxis in a microfluidic channel: mesoscopic and stochastic models

Chemotaxis is the response of a particle to a gradient in the chemical composition of the environment. While it was originally observed for biological organisms, it is of great interest in the context of synthetic active particles such as nanomotors. Experimental demonstration of chemotaxis for chemically-powered colloidal nanomotor was reported in the literature in the context of chemo-attraction in a still fluid or in a microfluidic channel where the gradient is sustained by a specific inlet geometry. In this work, we use mesoscopic particle-based simulations of the colloid and solvent to demonstrate chemotaxis in a microfluidic channel. On the basis of this particle-based model, we evaluate the chemical concentration profiles in the presence of passive or chemically active colloids, compute the chemotactic force acting upon them and propose a stochastic model that rationalises our findings on colloidal chemotaxis. Our model is also able to explain the results of an earlier simulation work that uses a simpler geometry and to extend its interpretation.


Introduction
Nanomotors are nano-to micro-meter sized machines that use a local (internal or found in the close environment) source of energy to move or perform work. The catalytic conversion of a chemical fuel on catalytically coated colloids was used with success for different types of motors (metallic rods [1], polystyrene [2] and silica [3] Janus particles, dimers [4]). Chemically powered nanomotors represent very promising devices for the execution of tasks such as sensing or cargo delivery in nano-to micrometer scaled environments [5][6][7]. Understanding the response of nanomotors to chemical concentration gradients, chemotaxis, is critical to engineer these possible applications.
The experimental characterisation of chemotactic motion can proceed either via chemical sources that supply the environment with fuel, at specific "target" locations, or via a flow that allows the sustainment of the gradient. The first idea was used by Hong et al to demonstrate the occurrence of chemotaxis for rod nanomotors [8]. The second strategy was used by Baraban et al in Ref. [9] where the authors studied experimentally the chemotactic motion of chemically powered nanomotors in a microfluidic channel and observed a lateral deviation of the nanomotors when the fuel is input asymmetrically at the inlet of the cell. The successful modeling and understanding of colloidal chemotaxis is also relevant for biological systems, with the traditional example being the chemotactic behaviour of bacteria [10]. The focus on biological chemotaxis has been put forward by Sengupta et al who also used a microfluidic channel to observe the enhanced migration of enzymes across a microfluidic channel in the presence of a * ORCID:0000-0003-4488-1615 † Corresponding author ‡ ORCID:0000-0002-6640-6463 substrate gradient [11]. More recently, the sorting of enzymes on the basis of their chemotactic response has been demonstrated in Ref. [12]. Numerical simulations of chemotaxis can either assume an expression for a "chemotactic force" or reproduce the chemotactic process itself, using a direct simulation of the solvent. This article starts with the latter approach, using Molecular Dynamics simulation with an explicit solvent in which several chemical species are represented. We then formulate a continuous picture for the solvent concentration field and use it to build a stochastic model that captures the chemotactic behaviour and its connection to our simulation parameters.
The literature on particle-based modeling of chemotaxis for nanomotors is rather scarce. Chen et al [13] considered a simplified system in which a constant chemical gradient is imposed by the boundaries of the simulation cell. In this work, we seek to imitate the experimental setup of Baraban et al [9] and to improve the theoretical understanding of this experimentally relevant configuration for chemically powered nanomotors, albeit using a simpler dimer-type nanomotor instead of the Janus and tubular microjet motors in the experiment. Our computational experiments builds gradually by starting with a non-catalytic spherical colloid, then adding a catalytic property to the spherical colloid and finally using the dimer nanomotor. We come back to the results of Chen et al [13] and extend their conclusions.
In section 2, we review the mesoscopic simulation model. The simulational implementation of the experiment of Ref. [9] is laid out in section 3. The continuous representation of the fluid's diffusion profile, including in the situation where chemical reactions occur on the surface of the colloid, is given in section 4. There, we also compute the chemical gradient induced forces on the colloids. The stochastic model for the colloids is presented in section 5. We give the results for both modeling strategies in section 6 and conclude in section 7. In appendix A, we provide full information on how to access the numerical code to reproduce our findings.

Simulation model
The solvent consists of point particles with a mass mi, a position ri and velocity vi. There is no force acting between solvent particles, their interaction is instead modeled via cell-wise collisions at fixed time intervals τ using the Multiparticle Collision Dynamics (MPCD) collision rule introduced by Malevanets and Kapral [14,15]. MPCD has been used successfully to investigate the dynamics of colloids in microfluidic channels by Prohm et al, in the context of inertial focusing [16], or by Nikoubashman et al to study the flow of colloids in the presence of obstacles [17], for instance.
The evolution of fluid particles, in the presence of forces due to colloids or to the flow-inducing field, is resolved numerically using the velocity Verlet algorithm [18,19] ri(t + dt) = ri(t) where fi is the total of pair forces on particle i and g is the external acceleration. The timestep dt for MD is a fraction of the one for MPCD dt = τ NMD , where NMD is the number of MD steps between successive MPCD collisions. At fixed time intervals τ , the fluid particles are sorted in a lattice of cubic cells of side a, whose origin is shifted randomly to ensure Galilean invariance [20], and their velocities are collided cell-wise according to where the prime denotes the post-collision value, ξ is the cell containing the particle i, Ω ξ a rotation operator of angle Θ in R 3 around a randomly chosen axis and v ξ is the centre-of-mass velocity of the cell. The transport properties of a MPCD fluid can be computed analytically [21,22, and references therein].
The simulation parameters are given in table 1. One technique to obtain a Poiseuille flow is to impose a constant acceleration field in the simulation cell [23], in combination with periodic boundary conditions in the direction of the flow (x) and stick boundary conditions in the direction transverse to the flow (z). In Refs. [19,24], this technique was applied to MPCD fluids. In the y direction we use specular boundary conditions for the fluid, so that the flow velocity profile vx(z) bears no dependence on y. We implement stick boundary conditions for the flow in the z direction with a combination of bounce-back collisions and ghost particles at the boundary MPCD cells [25]. The ghost particles also set the temperature at the wall. In this work the walls also provide a sufficient thermostatting action to compensate for the flow-induced heating and no bulk thermostatting is used. We have verified that the temperature of the fluid remains stable after a transient period at a value that is about 3% in excess of the temperature set at the walls.
The colloids evolve according to the velocity Verlet algorithm, similarly to the solvent particles, but do not participate in the collisions and are not subject to the acceleration field g. In the case of dimer nanomotors, the distance between the two spheres is held constant using the RATTLE algorithm [26].
The coupling with solvent particles is done via the shifted and truncated Lennard-Jones 12-6 potential, of the form for rij ≤ σij × 2 1/6 0 else (5) where ij denotes the strength of the potential and σij the radius of the colloid. i and j represent the species of the solvent and colloidal particles and rij is the distance between them. The variation of ij , depending on what type of solvent and colloid interact, leads to a net force in the presence of chemical concentration gradients. This will be made explicit in section 4. ij for interactions with the solvent of type A is set to 1 and defines the energy scale. In the following, all quantities are expressed in simulation units with length a of the unit cell, mass m of the solvent particles, energy A and time a 2 m/ A. As the interaction between all colloids and solvent particles of species A is equal, the first subscript to κ,A is dropped in the following.
Confinement of the colloids in the z direction is obtained by a purely repulsive Lennard-Jones 9-3 potential, of the form where z is here the distance to the closest horizontal wall. The same potential is applied in the y direction to avoid hypothetical crossings of the lateral wall during a simulation. MPCD fluids are well suited to describe chemical reactions influenced by catalytic surfaces [27] or in the bulk [28]. At the surface of a catalytic colloid C, the reaction converts fluid particles of type A (the fuel) to fluid particles of type B (the product). The reaction occurs when the fluid particles crosses the interaction region of the colloid and is executed with a probability p ∈ [0, 1], when the fluid particles exits the interaction region to avoid any discontinuity in the energy [27]. All simulations were performed with the open-source RM-PCDMD software [29,30].

Cell design
The inspiration for the design of the microfluidic cell comes from the experimental work on the chemotaxis of nanomotors by Baraban et al [9]. There, a channel is fed with fluid at a fixed flow rate via three inlets. As only an inlet contains a chemical species of interest, i.e. the fuel for the nanomotors, the lateral distribution of species is inhomogeneous. The effect of the resulting gradient is a chemotactic behaviour that is observed by monitoring the deviation angle of the nanomotors with respect to a straight line motion.
To reproduce the experimental features in a simulation, we have setup a thin channel with a forced flow in the program chemotactic cell of the RMPCDMD software that is illustrated in Fig. 1. A simulation snapshot in Fig. 2 shows the Poiseuille flow, the concentration field for the fluid species A and the trajectory of a dimer nanomotor. To avoid reproducing the inlet channels at a large computational cost, we assign the solvent species to be either A (the fuel) or F (a neutral fluid species) arbitrarily, depending on the lateral position of the particles. The total number of particles in the simulation is constant and particles are not created nor destroyed but "recycled" when they re-enter the simulation cell due to the periodic nature of the x boundary. The region of the cell where the forced attribution of species is applied is called the buffer and is located in the region 0 < x < L buffer . In the following of the paper, the trajectories from the mesoscopic simulations are shifted in x by −L buffer to match the coordinates system of the stochastic model. Figure 1: Schematic description of the simulation cell for chemotaxis in a flow. There are two inlets on the left. In the buffer (leftmost region), solvent particles are set to species A for 0 < y < Ly/2 and F else. Lennard-Jones 9-3 potentials confine the colloids close to z = Lz/2. The colloid are initially placed in the F inlet and constrained to a fixed y and z track. This constraint is released when x > σ. The confining plates are made of transparent blue. The dimer is made of a catalytic (red) bead and of a non-catalytic (blue) bead and its centre-of-mass trajectory is shown as a white line. Here, the nanomotor is attracted toward the front of the cell (lower values of y) and displays a corresponding reorientation. The instantaneous concentration field of A solvent particles around z = Lz/2 is shown in pseudocolor.
Although the simulations are three dimensional, the motion of the colloids is limited around the centre of the cell in the zdirection by the confining walls.
The flow between two plates is of the Poiseuille type with maximal velocity v flow and average velocity vav = 2/3 v flow . We report the characteristic numbers of the flow in Table 2. For the work of Baraban et al [9], we use the values found in the article for the flow rate (140µL per hour), the width (200µm per in-let, there are three inlets) and temperature (300K). The height of the channel was confirmed by email to be 30µm. For the properties of water, we use reference data from the NIST [31, retrieved January 11, 2017] (density ρwater = 996.56 kg/m 3 , and viscosity ηwater = 8.54 10 −4 Pa·s) and take the diffusion coefficient Dexp = 2 10 −9 m 2 /s that is a reasonable approximation for both water [32] and hydrogen peroxide [33]. For the MPCD fluid properties, we refer to the review by Kapral [22]. The viscosity is and the self-diffusion coefficient is We have verified that the velocity profile obeys the theoretical prediction [24] vx We have aimed for a similar fluid regime with respect to the experiment, that is a high Péclet number (Pe) flow, as our estimate for the concentration profile relies on Pe 1 (see section 4 and Ref. [34]). A value as high as for the experiment could not be obtained, but the important feature is that the transport by the flow dominates the one by diffusion in the x direction. The Reynolds number Re should remain moderate for the laminar regime to hold.
The Péclet, Reynolds, and Mach number for the flow are computed as Pe = Lzvav/D , m is the speed of sound [16]. The maximum Mach number is obtained for the maximum velocity of the flow vmax = ρ g L 2 z 8η ≈ 0.095 and is Ma ≈ 0.13. Besides the geometry of the cell, the simulation protocol differs slightly from its experimental counterpart. The colloid moves on a track at fixed y = Ly 2 + y shift and z = Lz 2 . This restriction is lifted when the x position of the colloid (centre of mass) has passed L buffer plus its own radius. This allows for a systematic comparison of the chemotactic drift across repeated runs without suffering from possible disturbance from the resetting of particles in the inlets. The shift y shift in the y direction ensures that solvent particles in the interaction range of the colloid are only of species F and not influenced to a change of species in the buffer region.

Density profiles and surface interaction
In this section, we compute the stationary concentration field cα(r) for the different chemical species in the cell and the resulting chemotactic force on the colloids.
We approximate the evolution of cα(r) by a 1D diffusion equation where the spatial direction of the flow is proportional to time, i.e. x = v flow t, and the diffusion process acts on the transverse direction y. This approximation has been tested experimentally in Ref. [34] and is only valid close to z = Lz/2. Close to the boundaries z = 0 and z = Lz of the cell, a different transport regime is taking place. Here, we neglect the variations in the z direction for the computation of the concentration cα(r) in the absence of catalytic particle.
Fixing z = Lz/2 and identifying the time and the x coordinate results in The separate inlets for the A and F chemical species translate into the initial condition cA(0, y) = c0Θ(Ly/2 − y), cB(0, y) = 0 and cF (0, y) = c0Θ(y − Ly/2) for x = 0 (equivalently t = 0). The solution to Eq. (11) is cB(x, y) = 0 (14) where erfc is the complementary error function. In the absence of a catalytically coated colloid, the value of cB remains zero at all times. The average number density is a constant, in the low Mach number conditions here.
In the remainder of this section, we use coordinates centered around the colloid. The spherical coordinates are defined by the following relations    x = r cos ϕ sin θ y = r cos θ z = r sin ϕ sin θ (16) and are represented in Fig. 3. The presence of a catalytic colloid is taken into account in Eq. (11) by a radiation boundary condition (RBC) on the surface of the colloid, at radius R. The boundary condition is applied at the limit of the interaction region, i.e. R = 2 1/6 σ, where the continuum diffusive picture breaks down and where the chemical reactions are triggered. The RBC identifies the flux of chemicals with the consumption of the catalytic reaction, at the surface of the colloid. It was developed by Collins and Kimball [35] and rederived later by other authors [36,37]. The RBC at radius R in the absence of external gradient, for the species A, is where kD = 4πRD is the diffusion-limited rate constant and k0 = pR 2 8πkBT /m is the chemical rate [38]. p is the reaction probability defined in section 2 and m is the mass of a fluid particle.
x y z θ φ r Figure 3: The spherical coordinates centered on the colloid.

Surface interaction
Following the work of Rückner and Kapral [27], we sum the interaction potential Vκ,α between a colloid of species κ and the fluid particles of species α over the interfacial region to obtain the potential energy. Differentiation with respect to the colloid's coordinate gives the force on the colloid: In all the subsequent derivations, we rely on the solution of the diffusion equation in the bulk of the fluid. Within the interfacial region, the interaction with the colloid modifies the concentration cα by the Boltzmann weight [27]. Equation (18) where β = (kBT ) −1 and 1r is the unit vector pointing in the radial direction. Integrating by parts yields where we have defined Within this framework, it is sufficient to compute surface integrals over the colloids to compute the force acting on them. The use of Eq. (20) in the literature is however limited to situations where there is no external chemical gradient.

Single passive colloid
For a single passive colloid of type N , we use the solution (12)- (14) together with Eq. (20).
Locally around the colloid, we use a first-order expansion of cα in the lateral direction.
where λ = ∂ycA(x, y) The explicit dependence of FN,y and λ on the coordinates is omitted for brevity.

Single active colloid
To take into account the catalytic activity of a colloid, the diffusion equation (11) is solved with the radiation boundary condition (RBC). Equation (17) expresses the flux that originates from the chemical activity of the colloid, that is the only contribution to the flux in the situation where the RBC was derived.
The total radial flux is the sum of the reaction-induced flux from Eq. (17) and of the diffusive flux 1 4πR 2 D 1r · ∇cA = RkDλ cos θ.
We make the following ansatz for the concentration field cA: Eq. (26) is a solution of the diffusion equation and matches the boundary condition (25). The first three terms come from a truncation of the Legendre polynomial expansion of a diffusion profile with a spherical catalytic sink. The last term is needed to reflect the presence the external gradient. We obtain the coefficients by inserting Eq. (26) in Eq. (25): where we have used for c0 the solution (12). The value of λ is obtained by differentiating Eq. (12) with respect to y.
It is important to mention that (i) in the absence of a gradient, the coefficients in Eq. (27) lead to the standard solution of a spherical sink and (ii) in the absence of chemical activity, the solution (12) is recovered to linear order.
To obtain the solution for cB, we observe that the reactive flux of A particles absorbed at the surface of the sphere equals the one of B particles, except that the currents flow in opposite directions. The solution to the diffusion equation with these opposite flows on the surface is To compute the force in Eq. (20), we set ΛC,F = ΛC,A. Using Eq. (15), we find that the chemotactic force on the active colloid C is 1 The flux is multiplied by the surface of the sphere 4πR 2 as this is how the RBC is expressed in the literature.

Dimer nanomotor
A dimer nanomotor is made of two spheres linked rigidly, with a distance d between their centre of masses. One (C) is catalytic and acts as a sink for A-species solvent particles and a source of B-species solvent particles, creating locally a gradient centered around C.
For the dimer nanomotor, the concentration fields depend on the presence of the catalytic sphere of the motor in the same way as for the single active sphere and we reuse Eqs. (27) and (29). The force must be evaluated also on the non-catalytic sphere N where the spherical symmetry does not hold. We evaluate this expression numerically as where the prime denotes the spherical coordinates around the N sphere of the dimer, illustrated in Fig. 4.
The forces on C and N are summed to obtain the centre-ofmass (com) force: The torque on the dimer is The motion of the dimer is studied via its centre-of-mass position in the x-y plane and its inclination φ with respect to the x axis.

Langevin dynamics for the colloids
For the spheres, the evolution of x and y (here, coordinates in the laboratory frame of reference) is given by the overdamped Langevin equations [39].
where v flow is the average velocity due to the flow at z = Lz/2, ξx and ξy are normal distributed white noise, D is the diffusion coefficient of the colloid and γ its friction. For the single spheres, we use γ = 4πησ, the slip hydrodynamics friction, and D = kBT /γ.
For the dimer, we use the following Langevin equations for the centre of mass and orientation where the projected forces F and F ⊥ are and T is the torque defined in Eq. (32). The rotation operations in Eq. (35) and (37)  is an anisotropic colloid with axial symmetry. Accordingly, we compute separately the diffusion coefficient parallel ( ) and transverse (⊥) to its axis as where ζ is or ⊥, following Ref. [39]. The projected velocities are defined by withû the unit vector in the direction joining the two spheres. The angle φ measures the deviation ofû in the x-y plane with respect to the unit vector 1x. The results of these simulations is D = 2.0 10 −3 and D ⊥ = 1.5 10 −3 . The rotational motion is characterised in the same fashion, using the velocity autocorrelation of the angle φ: φ (τ )φ(0) with the result Dr = 1.4 10 −4 . We determine the friction coefficient using the fluctuation-dissipation relation The rotational time D −1 r ≈ 7100 is such that in the chemotaxis simulations the dominant angular behaviour reflects the rotational drift. The contribution of rotational diffusion is nevertheless visible in the spread of angles in Fig. 10. A typical trajectory for the passive sphere is represented in the x-y plane in Fig. 6. The main component is a displacement to the right under the influence of the flow, to which thermal fluctuations are superimposed. For this trajectory, N,F = A so there is no chemotactic behaviour.
Upon changing N,F , a lateral force is exerted on the colloid due to the combination of the asymmetry of cA and cF with respect to the colloid and of the difference in surface interaction between the colloid and the A and F solvent species. This situation is one of passive diffusiophoresis.
The results for passive spheres are summarised in Fig. 7     Conversely, for N,F > A, the chemotactic drift is negative. This is shown for N,F = 2 and N,F = 4 in Fig. 7. Figure 7 reports the results for both the mesoscopic model and the stochastic model. Given the fluctuating nature of the colloid dynamics, we have repeated the simulations 16 times and denote by a shaded area one standard deviation below and above the average results. We observe that the average trajectories are similar and also that the spread of trajectories is similar.
This first result on passive chemotaxis provides here the simplest context for the experimental setup and allows a possible calibration of the quantity ΛN,A − ΛN,F without any dependence on reaction rates.

Active sphere
We now turn to the chemotactic behaviour of an active sphere. The injection setup is the same as for the the passive sphere but now C,F = A and the mechanism that creates a systematic lateral motion for the passive sphere, that is proportional to ΛC,A − ΛC,F , is effectively zero.
What occurs, instead, is that the local gradient in cA generates an asymmetric distribution also for cB that, together with a nonzero value for ΛC,B −ΛC,A, creates an original combination of passive and active diffusiophoresis. Using Eq. (29) and observing that c2 > 0 (as λ < 0 everywhere in the simulation), we expect a downward chemotaxis for ΛC,B − ΛC,A > 0 (i.e. C,B < A) and an upward chemotaxis for ΛC,B − ΛC,A < 0 (i.e. C,B > A). We explore the effect of changing C,B in Fig. 8 that confirms this direction for the chemotaxis of the active colloid. As for the passive sphere, the comparison between the mesoscopic and the stochastic models is positive.
Even though the situation of the active sphere is simpler than the experiments on nanomotors of Ref. [9], it already contains a complex ingredient: the chemotactic force is entirely caused by the self-generated concentration field around the colloid. Thanks to the derivations in section 4, we understand how the asymmetry of the imposed concentration fields of A and F lead to an asymmetry in the concentration of B that gives rise to chemotaxis.

Dimer nanomotor
For the dimer nanomotor, we track in the simulations the centreof-mass position rcom and the orientation φ. The release from the injection track occurs here as soon as both spheres have exited the buffer region and the initial orientation is φ = 0.
We have chosen, for the simulation parameters, that C,α = N,α, for all solvent species α. This is at variance with simulations of the prototypical dimer nanomotor [4,13,27,40] where C,A = C,B = A = N,B . Preliminary tests showed that in this situation the trajectories are pathological in the sense that there was no "gentle" deviation of the orientation φ and that no systematic tendency could be found. For nanomotors in a bulk environment with no a priori gradient in the chemical concentrations, this choice would be irrelevant.
Both the lateral deviation in the x-y trajectories and the angular deviation for φ are reported in Figs. 9 and 10 and are consistent between the mesoscopic simulations and the stochastic model. This confirms the adequateness of the combined forcetorque evolution given by Eqs. (35)-(36) as both the direction in y and in φ are reproduced. The origin of the lateral deviation can be understood as coming not only from the redirection of the selfpropelling force, that is oriented along the dimer axis, but also from the net lateral force on the centre of mass whose direction we can infer from Eqs (27), (29) and (31). To our knowledge, this is the first time that the continuum force computation of Ref. [27] is extended to compute the lateral force on a nanomotor and the resulting torque on the dimer nanomotor.

Comparison to earlier simulation work
Aside from the geometric considerations that are specific to the microchannel, the stochastic model that we have designed can accommodate other situations. We use it to test, in a qualitative sense, the observations of Ref. [13]: we re-use the diffusion coefficients that we have obtained with our cell geometry and fluid parameters but change the type of gradient and the interaction type of the motor to match those of the earlier paper, which we call the "constant gradient" model. In Ref. 2. The average trajectory of the nanomotor is toward y = Ly (see Fig. 7(a) of Ref. [13]).
We will show here that these results depend strongly on the choice of parameters for the interaction between the colloid surfaces and the solvent via three observations: the distribution of trajectories y(t), the distribution of angles P (θ), and the distribution of position P (y). We have performed stochastic simulations with κ,A = κ,F = C,α = 1 (for all solvent species α and colloid species κ) and N,B = 0.25, Ly = 30, Lz = 15. A harmonic wall with spring constant k wall = 10 prevents the exit of the colloid and is turned 2 For clarity, we reuse the species labels and orientation of the gradient used in the present work. This should be kept in mind when comparing with the work of Chen et al where the gradient is along x and the species are labelled F for the fuel (here, A), S for the inert fluid (here, F ) and P for the reaction product (here, B). on for y < 2 or y > Ly − 2. The nanomotors were started in the x direction, as in Ref. [13], to avoid an a priori bias. The results are displayed in Fig. 11. The distribution of y(t) trajectories in the upper panel shows a large spread, from which we cannot conclude of a dominant chemotactic behaviour. Examining P (y), we see however that there is an average accumulation of particles close to y = 0. The angular distribution P (θ) does show a bias toward y = 0 (equivalently θ = π), as in Ref. [13].
Until here, there was no chemically-induced force on the C sphere. To assess the importance of this choice, we perform another set of simulations C,B = N,B that are shown in Fig. 12. This parametrisation is the one used for our main results of sections 6.1, 6.2, and 6.3. From Eqs. (29)-(32), we know that the torque on the nanomotor will be influenced by this choice. This is reflected in the angular distribution P (θ) in Fig. 12 that now displays a strong orientation toward θ = 0 (i.e. the nanomotor is oriented toward y = Ly). As a result, there is no competition between a downward-facing nanomotor and the greater velocity for upward-facing nanomotor orientation that is observed in Ref. [13] and the distribution of y(t) trajectories is narrow, showing a strong chemotactic behaviour. This is confirmed by the distribution of positions P (y).
To conclude this comparison, we have observed competing effects of orientation and propulsion, leading to different chemotactic behaviours. The distribution of angles, for instance, depends stronly on the parameters chosen for the colloid-fluid interaction: For the interaction choice of Chen et al, C,A = C,F = C,B , we also find an average orientation toward y = 0, but for the other choice that we used the average orientation is toward y = Ly. The overall chemotactic behaviour will however also depend on the diffusive and propulsive properties of the dimer.
Although further research on the experimental characterisation of the surface properties of the nanomotors is needed, our stochastic model provides an interesting tool to relate these properties to the effective chemotactic behaviour of nanomotors.

Conclusions
We have proposed a particle-based simulation setup, based on Multiparticle Collision Dynamics and Molecular Dynamics, to study the chemotactic motion of passive and active colloids. The concentration field that drives this motion is sustained by the flowing input of the cell, via two inlets, as is done in the experiment of Ref. [9]. We then constructed an approximate solution for the chemical concentration profile in the microfluidic channel in the presence of an active colloid and built a stochastic model with a microscopic expression for the systematic force. This model provides a sound basis for the microscopic origin of chemotaxis, i.e. it explains why the colloids move towards higher or lower values of the coordinate y. Upon extending this model to a two-sphere dimer nanomotor, we also gain an understanding of why the nanomotor changes its orientation via a systematic torque. We have thus improved on the formal understanding of the chemotactic motion in the microfluidic channel. Our stochastic model, adapted to the setup with a constant concentration gradient by Chen et al reproduces the observation that the nanomotor tends to orient opposite to the gradient. While Chen et al observed a net positive chemotaxis, this is not the case in the present work. The reason is that the interplay between orientation and propulsion depends on the geometry of the motor and on the choice of parameters. This reasoning is supported by changing the surface interaction also for the C bead, resulting in a cooperation of orientation and propulsion leading to positive chemotaxis. The specific change in C,B , that is not used in the simulation literature but introduced here in subsection 6.3, allowed us to probe a regime of well-defined positive chemotaxis and demonstrates that the class of simulation models introduced by Rückner and Kapral [27] possesses a very rich phenomenology.
It is interesting to note that our stochastic model requires only the input of the surface interaction parameters Λκ,α and of the diffusion coefficients of the dimer. In principle, it could be extended to other forms of the interaction potential or rely on a characterisation of Λκ,α that does not reveal the full functional form of this potential, and also to other motor geometries. While we have chosen ΛN,α = ΛC,α for all solvent species α for this first investigation of the model, there is room for possible qualitative changes with other choices, as we have witnessed for the constant gradient configuration.
The overall speed of the motor, as it was the case in earlier simulation studies for dimer nanomotors, only depends on ΛN,A − ΛN,B [27]. The torque, on the other hand, depends also on ΛC,A − ΛC,B. A difference in these quantities, and thus in the surface properties of both sides of the nanomotor, could be revealed in a controlled manner in experiments and this reinforces the importance of the microfluidic channel configuration for chemotactic studies.
The present work can be extended to other types of motors, notably the Janus nanomotors that are used in Ref. [9] and for which colloidal assemblies have already been used in mesoscopic simulations [41]. This line of research is promising to test in silico the behaviour of different motor geometries. Using models suitable for the chemo-mechanics of enzymes, at a mesoscopic level [42][43][44], could provide very fruitful advances for understanding the recent works on enzyme chemotaxis [11,12], especially given the fact that multiple inlets microfluidic devices originate from studies on bacterial chemotaxis [45] and have been used for the enzyme studies [11,12].