1. Introduction

In 1970 Arthur Ashkin first demonstrated how to manipulate and confine microscopic particles suspended in water through radiation pressure [1]. Following the first demonstration of optical trapping, Ashkin and collaborators developed the single-beam gradient trap, today known as optical tweezers (OT) [2,3]. The basic principles of OT utilise the fact that light carries momentum which can be harvested to manipulate microscopic particles in solution. In its conventional and simplest set-up, OT focus a collimated Gaussian laser beam to a diffraction-limited spot where it can trap microparticles. Soon after the first demonstration of OT, Ashkin et al. employed them to manipulate biological particles like bacteria and erythrocytes without causing damage [4,5].

In humans, erythrocytes, or red blood cells (RBCs), are anucleated cells responsible for the oxygen delivery to tissues and organs. Mature and healthy RBCs have a biconcave disk shape that minimizes the membrane bending energy. Typically, RBCs have diameter of 6-8 $\mu \rm {m}$, a peripheral thickest portion of 2-3 $\mu \rm {m}$, and a central dimple 0.8-2 $\mu \rm {m}$ thick [6]. The excess surface area and membrane elasticity render the cell elastic and permits the RBC to pass through the microvasculature by deforming [7]. Alterations in the RBCs’ membrane elasticity are implicated in severe disfunctions of the microcirculation (e.g., capillaries can be entirely clogged, triggering tissue necrosis or organ damage and failure) [8]. The RBCs’ alteration can be genetically inherited [9], a consequence of a pathogen infection [10], a metabolic disorder [11], or due to radiation treatment [12]. Very recently, it has also been correlated to SARS-Cov2 infection [13].

In the last decades, OT have been widely applied in RBC research to investigate biochemical and biophysical properties of both healthy and unhealthy RBC via single- or multi-beam OT [14]. In these studies, researchers have adopted two main approaches to trapping: the indirect trapping, where handles as silica or polystyrene microspheres are used to manipulate the RBC [15], and the direct trapping, where the light beam directly traps the RBC [5]. Studies on the RBC optical trapping have demonstrated that high laser power (>500 mW) induces deoxygenation of the trapping site, cell membrane ionization, inactivation and cells expulsion from the trap by the radiometric force. Moreover, the undesirable temperature rise due to the incident laser light (of the order 1.4 C$^{\circ }$/100 mW [16]) which affect the Brownian motion adding noise to the system and thus setting a limit of resolution. Also, under high laser power the shape change induced by the radiation pressure could cause strong temperature gradient across the cell membrane which could lead to membrane rupture [17]. Therefore, regardless from the mechanism of trapping, the nature of biological samples makes them particularly susceptible to photodamage. To minimize this, infrared light in the second biological window (wavelength around 1064 nm) and minimizing optical power is generally preferred for experiments [18].

As the cell is significantly larger than the incident wavelength, the geometrical optics approximation (GO) models properly the beam cell interaction [19–21]. GO assumes that the beam can be discretized in a series of rays that carry a fraction of the total momentum and by calculating and summing up the scattering of all rays, it is possible to compute the total force applied. However, even though GO simplifies considerably the theoretical treatment compared to a full wave optical approach [22], an accurate calculation requires consideration of a large number of rays, with an associated high cost in computational time. Simulating the Brownian dynamics requires repetition of the force calculation at each time step sequentially, which becomes prohibitively slow if the force calculation is not optimized [23]. The fact that the calculation is sequential and that the shape is complex prevents the use of conventional approaches to speed up the calculation (e.g., parallelization and interpolation based approaches).

Machine learning, and in particular neural networks (NN), are emerging across a variety of research fields as a powerful technique to solve challenging problems. Backed by their ability to learn from previous examples in order to make new predictions, NN are contributing to biology [24], food sensing control [25], and even to containment of epidemics [26]. In fact, NN have recently been demonstrated as an useful technique to increase both the speed [27] and the accuracy [28] of optical forces calculations when compared to GO, allowing the study of more complex systems through Brownian dynamics simulations. While these previous works consider spheres [27] and ellipsoids [28], there is no evident reason to remain constrained to these relatively simple shapes. Indeed, the computation time saving that could be achieved by using NN for force calculation of particles with more complex shapes makes this a particularly attractive application.

In this work, we train a NN to enhance the speed and accuracy of the optical force calculation for RBC. This permits a numerical exploration of the Brownian dynamics of a RBC, potentially allowing to study in a more complete manner different trapping configurations. More efficient trapping configurations employ less laser power and therefore reduce the risk of photo damaging the trapped cells.

2. Methods

2.1 Model and geometrical optics calculations

In our model we consider a Gaussian beam propagating along the opposite direction of the force of gravity ($+z$ direction). The wavelength ($1.064 \, \mu \rm {m}$) is selected to match the vast majority of experiments and in agreement with previous works on trapping RBCs [11,17]. The OT parameters are the ones of a typical OT experiment (beam power $5 \, \rm {m}W$, numerical aperture $1.3$).

The RBC is assumed to be in its healthy biconcave disk conformation, and the parameters describing the shape of the RBC are those reported by Evans et al. [6]. A radius ($r$) of $3.91 \, \mu \rm {m}$, a central dimple with a thickness ($t_{min}$) of $0.81 \, \mu \rm {m}$, and a thickest portion, located at $2.76 \, \mu \rm {m}$ from the axis of symmetry, with a thickness ($t_{max}$) of $2.52 \, \mu \rm {m}$. According to the Evans-Fung model, the thickness ($Z$) of a section of the RBC reads:

As RBCs are significantly larger than the wavelength of the incident light, the optical forces acting on them can be calculated with GO. We perform this calculation with the specialized software OTGO [30]. For biological samples, such as RBCs, that have a low refractive index contrast with the typical suspending medium, the fraction of power that is reflected after a scattering event is very low ($< 0.001$) [31], therefore in our ray tracing calculations, only the first two scattering (refraction) events are considered.

2.2 Diffusion tensor

The erratic motion of a particle trapped in liquid in an OT set up is influenced by the fluid’s resistance, by the thermal noise, and by the external deterministic forces exerted by the OT [22,32]. For non-spherical objects, a single scalar diffusion coefficient is not enough to describe the statistics of the random motion. It is necessary to use a 6 x 6 diffusion tensor ($\mathbf {D}$), which depends on the particle shape and orientation [22]:

Although an analytical expression for ($\mathbf {D}$) exists for simple shapes like spheres, ellipsoids or cylinders [33], the RBC morphology is more complex and requires numerical methods for its determination. Here, we used the bead model technique developed by De La Torre et al., exploiting the widely used software winHYDRO++ [34,35]. In the bead model, a series of spheres are used to approximate the size and the total volume of the RBC. From the bead model, winHYDRO++ calculates the 6x6 tensor ($\boldsymbol \Xi$) encoding the hydrodynamic resistance of the non-spherical particle. We then obtain the diffusion tensor $\mathbf {D}$ via the generalised Einstein relationship [36]:

In the present study, the bead model is constructed in a strict sense, filling the volume of the RBC considering only spheres of equal sizes. In this case, the centre of diffusion of the particle coincides with the centre of mass of the particle, and the numerical output for the diffusion tensor reads:

2.3 Particle dynamics simulation

The simulation of the dynamics of the RBC is based on the works of M. X. Fernandes et al. [32], and described in Ref. [22]. Two reference frames are defined: a particle reference frame $\Sigma _{\rm {p}}$, which has an origin that coincides with the particle’s centre of mass (CM) and the centre of diffusion (CD), and a laboratory reference frame $\Sigma _l$ that is centred at (0,0,0) and which axes are oriented along $\mathbf {\hat {x}}$, $\mathbf {\hat {y}}$ and $\mathbf {\hat {z}}$, Fig. 1(a). At time $t$, the RBC’s CD is located at $\mathbf {r_{\rm {CD}}}(t)=[x_{\rm {CD}} (t),y_{\rm {CD}} (t),z_{\rm {CD}} (t)]$. The cell orientation can be described by the angles $\alpha _1 (t)$, $\beta _1 (t)$ and $\gamma _1 (t)$ defined with respect to the particle unit vector $\mathbf {\hat {x}_{\rm {p}}} (t) =[\hat {x}_{\rm {p,x}} (t),\hat {x}_{\rm {p,y}} (t),\hat {x}_{\rm {p,z}} (t)]$, $\mathbf {\hat {y}_{\rm {p}}} (t) =[\hat {y}_{\rm {p,x}} (t),\hat {y}_{\rm {p,y}} (t),\hat {y}_{\rm {p,z}} (t)]$, $\mathbf {\hat {z}_{\rm {p}}} (t) =[\hat {z}_{\rm {p,x}} (t),\hat {z}_{\rm {p,y}} (t),\hat {z}_{\rm {p,z}} (t)]$. $\mathbf {D}$ is obtained in the particle reference frame, that is centred at $\mathbf {r_{\rm {CD}}}(t)$ and the axes are oriented along $\mathbf {\hat {x}_{\rm {p}}}$, $\mathbf {\hat {y}_{\rm {p}}}$ and $\mathbf {\hat {z}_{\rm {p}}}$

 

Fig. 1. a) Definition of particle (yellow) and laboratory (white) reference frames and rotation angles ($\alpha$, $\beta$, $\gamma$) of the cell around the laboratory reference frame; b) schematic depiction of the neural network. The input layer contains six neurons describing the cell position and orientation, and the output layer has six neurons describing the components of force and torque acting on the cell. In between are seven hidden layers ($i=7$), each of them with 256 neurons ($j=256$). c-d) Density plots comparing the magnitude of the total force ($F_{\rm {tot}}^{\rm {NN}}$) and torque ($T_{\rm {tot}}^{\rm {NN}}$) predicted with NN with those calculated with the GO method ($F_{\rm {tot}}^{\rm {GO}}$) and ($T_{\rm {tot}}^{\rm {GO}}$). Regression lines are shown in red. e) Log-Log plot of the normalised root mean squared error (NRMSE) between $F_{\rm {tot}}^{\rm {NN}}$ and $F_{\rm {tot}}^{\rm {GO}}$, and $T_{\rm {tot}}^{\rm {NN}}$ and $T_{\rm {tot}}^{\rm {GO}}$ as a function of the number of rays used in the GO calculation. For each data point, the NN employed remains the same (trained with $4{\times }10^{2}$ rays).

Download Full Size | PDF

To simulate the free diffusion of an arbitrarily shaped particle from time $t$ to the time step $t+\Delta t$, initially one has to calculate the increment of the particle position and orientation in $\Sigma _{\rm {p}} (t)$:

Therefore, the finite difference equation to update the particle position in $\Sigma _l$ is:

Once the new particle position is calculated, one has to update the particle orientation from $\Sigma _{\rm {p}} (t)$ to $\Sigma _l (t)$, which is effectively a rotation of the particle unit vectors. This rotation, for small angles, is expressed in $\Sigma _{\rm {p}}$ by the rotation matrix:

Transforming this rotation matrix to $\Sigma _l$, we obtain the unit vectors representing the orientation of the particle at the end of the time step:

As the last step, the rotation matrix has to be updated:

However, in the current situation, we must also account for the optical forces ($\mathbf {F}$) and torques ($\mathbf {T}$) exerted by the optical trap on the centre of mass of the RBC. Therefore, taking into account $\mathbf {F}$ and $\mathbf {T}$, the increments of the particle orientation and position in $\Sigma _{\rm {p}}$ are:

The finite difference scheme is combined with the NN or with the GO code to calculate the optical forces and torques acting on the RBC and to simulate the Brownian dynamics of the optically trapped particle. We estimate the time step, $\Delta t$, for the Brownian motion simulation from the trap stiffness reported by Tognato et al. [21], and from the diffusion properties of a healthy RBC obtained here. The typical time scale on which the restoring force acts is given by $\tau _{\rm {OT}}=\frac {\gamma }{k}$, while the momentum relaxation time is given by $\tau _{\rm {m}}=\frac {m}{\gamma }$, where $\gamma$ is the particle friction coefficient, $m$ is the meass of the particle and $k$ is the trap constant. To assure numerical stability, $\Delta t$ must fall in between these two characteristic time scales ($\tau _{\rm {OT}}\gg \Delta t\gg \tau _{\rm {m}}$) [22]. From the diffusion tensor $\mathbf {D}$, one can extract the diffusion properties of the RBC along a specific direction ($\rm {D}_i$), then through the fluctuation-dissipation theorem one can obtain $\gamma _i$. For example, for the x-direction one obtains $\gamma _{\rm {x}}=\frac {k_B \rm {T}}{D_{\rm {x}}}\sim 5{\times }10^{-8} \frac {\rm {kg}}{\rm {s}}$. Therefore, considering a trap stiffness in this direction $k_{\rm {x}} \sim 1.6 \frac {\rm {pN}}{\mu \rm {m}}$, one obtains $\tau _{\rm {OT}} \sim 3{\times }10^{-2}\rm {s}$. On the other hand, given a mass of $\sim 1{\times }10^{-11} \rm {kg}$ for a typical healthy RBC, one obtains $\tau _{\rm {m}} \sim 4{\times }10^{-4}\rm {s}$. A similar estimation can be made for the other directions, and contemplating the magnitude of the other terms in $\mathbf {D}$, a time step $\Delta t=0.001\rm {s}$ is adequate to assure the numerical stability [22].

2.4 Neural network architecture and training

The neural network (NN) architecture is composed of one input layer with 6 neurons representing position and orientation of the cell ($x$, $y$, $z$, cos$(\theta )$, sin$(\theta )$, $\phi$), one output layer with 6 neurons representing the force and torque components ($F_{\rm {x}}$, $F_{\rm {y}}$, $F_{\rm {z}}$, $T_{\rm {x}}$, $T_{\rm {y}}$, $T_{\rm {z}}$), and 7 hidden layers in between with 256 neurons each, Fig. 1(b). While $\Sigma _{\rm {p}}$ encodes the orientation of the particle by using 3 angles, because of symmetry, the orientation of the RBC can be completely defined by the polar angle $\phi$ and the azimuthal angle $\theta$, as shown in Fig. 1(a).

The total data generated with GO calculations is composed of $6{\times }10^{6}$ points. We divide this data into three data sets: One for training ($80{\%}$ of the data), another one for validation ($10{\%}$ of the data), and a last one for testing the final network performance ($10{\%}$ of the data). The training data are generated via GO calculations made in OTGO [30]. The cell is placed in uniformly distributed positions in a cube of side $8 \mu \rm {m}$ centred at the origin of the Cartesian coordinates system (i.e. $-4 \mu \rm {m} \leq x \leq 4 \mu \rm {m}$, $-4 \mu \rm {m} \leq y \leq 4 \mu \rm {m}$ and $-4 \mu \rm {m} \leq z \leq 4 \mu \rm {m}$). Simultaneously, to account for the possible different orientation of the RBC within the trap, the cell is uniformly and randomly oriented in an interval for $-\pi \leq \theta \leq \pi$ and $0 \leq \phi \leq \frac {\pi }{2}$. The training data are generated for the simplest case of a single-beam OT with the geometrical focus centered at $(0,0,0)$ and a beam power of $5 \rm {mW}$.

The NN is trained in Python using Keras (version 2.2.4-tf) [38]. The training of the NN is divided into 5 different steps. The data pre-processing and the model definition, which are done only once, and the loading of the data, the training step, and the evaluation of the performance, that are carried out iteratively. The training data, generated as previously described, contains data in different units and scales. While the position scale is in the order of $\sim 1{\times }10^{-6} \rm {m}$, the forces are on the range of $\sim 1{\times }10^{-12} \rm {N}$, and the torques are around $\sim 1{\times }10^{-18} \rm {N}\cdot \rm {m}$. To achieve an efficient training of the NN, we need to apply a pre-processing step where the variables must be rescaled around unity and $\theta$, that ranges from $-\pi$ to $\pi$, is expressed in terms of sines and cosines to avoid inconsistencies around $2\pi$. Shuffling the data and dividing them into a validating and training set is the final step of the pre-processing. In our case, the training data set contains $4.8{\times }10^{6}$ points, $6{\times }10^{5}$ different points are used for the validation data set, and $6{\times }10^{5}$ points are reserved for the testing data set. In this work, we employ fully connected NNs where each neuron is activated by a sigmoidal function. Defining the model implies choosing the number of layers and the number of neurons per layer. Among the explored architectures, the one consisting of 7 hidden layers provides the best results (in terms of accuracy, training time, and speed).

The iterative part of the training starts by loading the training data and applying the training step where the NN weights are optimised to minimise the loss function. We use the mean squared error as the loss function and the Keras implementation of the Adam optimiser [38]. Once the training dataset is fully explored, the difference between the NN calculation and the validating dataset (defined as the mean square difference) is computed. The iterative step is repeated until this difference reaches its minimum value and we consider that the model is fully trained. The training of the NN is done in a GPU type NVIDIA GeForce RTX 2060 with 16 GB of memory. The processor of the computer is an Intel Core i7-10700, and it has 16 GB of RAM.

3. Results

3.1 Single beam optical tweezers

To evaluate the effectiveness of our approach, we start by testing the ability of the NN to predict the forces and torques acting on an RBC in a single beam OT (SBOT). We compare the NN predictions (trained with data generated using $4{\times }10^{2}$ rays) and the GO calculations considering 4 times more rays ($1.6{\times }10^{3}$ rays) at $1{\times }10^{5}$ random positions and orientations. The 2D density plots shown in Fig. 1(c and d) illustrate the agreement between the NN and GO in predicting the optical forces (regression coefficient $0.998$, $R^2=0.996$) and torques (regression coefficient $0.999$, $R^2=0.996$), respectively. We further demonstrate the accuracy of the NN by comparing our NN (trained with data generated with $4{\times }10^{2}$ rays) with the GO calculation (considering a greater number of rays). Figure 1(e) shows the normalised root mean squared error (NRMSE) between the predictions of the NN trained with $4{\times }10^{2}$ rays and the GO calculations with different numbers of rays (up to $5{\times }10^{3}$ rays). The NRMSE decreases as the number of rays increases. The forces and torques calculated with $5{\times }10^{3}$ rays result more similar to the NN output than to the forces obtained with a total of $4{\times }10^{2}$ rays, meaning that the NN is able to increase the accuracy of the force and torque prediction, even for an object with such a complex shape.

3.2 Double beam optical tweezers

Since the NN is trained for a SBOT, one may think it can only predict the optical forces and torques for a SBOT. However, the NN can be used multiple times to simulate multi-beam optical tweezers. In fact, the NN can predict the forces generated by a single beam on different locations on the cell, and the total force acting on the centre of mass of the cell may then be calculated as the vector sum of each contribution. The experimental implementation of a multi-beam OT setup presents greater challenges in the beam alignment, power balance and beam control compared to a single-beam configuration. However, recent advancements in the field of OT and beam shaping techniques have made it possible to realize the potential of multi-beam OT in experimental setups [22,39,40]. Here we consider a double-beam optical tweezers (DBOT) where the two beams’ geometric foci are positioned $5.06 \mu \rm {m}$ apart along the $x$-axis, similar to the experiments conducted by Agrawal et al. [11], Fig. 2(a). To the best of our knowledge, the cell configuration observed by Agrawal et al. is the only one observed experimentally when a double beam optical tweezer is employed for trapping. Indeed, optical torques and forces are responsible to maintain the positional and orientational equilibrium of the cell. In fact, for any displacements from the equilibrium configuration restoring torques/forces act on the cell pushing it back to the equilibrium position and orientation [21].

 

Fig. 2. a) Schematic depiction of an RBC trapped by a double-beam OT. b-c) Comparison between the GO calculation and the NN prediction for the b) torque-rotation curve for rotation around the x-axis and c) force-displacement curve along the x-direction. (d) Comparison of the probability distribution obtained with the GO calculation and with the NN prediction for a RBC in a DBOT. (e) Cell orientations in the numerical simulation for both GO and NN.

Download Full Size | PDF

Figure 2(b-c) shows $T_{\rm {x}} (\alpha )$ and $F_{\rm {x}} (x)$ calculated with GO and predicted with the NN for a cell in its folded configuration (i.e., cell major axis parallel to the optical axis) trapped in a DBOT. In both cases, the NN predictions (solid line) agree well with the GO method (dots), demonstrating the possibility to use the NN for multi-beam optical traps. We therefore conclude that this approach can be extended to predict forces and torques generated by a three- and four-beam OT, situations in which the GO calculation is considerably slower given the very large number of light rays required.

We now investigate the cell’s dynamics within a DBOT using both NN and GO to compute the optical forces. The simulation of the Brownian dynamics follows the strategy explained in the Methods section (Particle dynamics simulation) where now, the force and torque considered is the sum of the contributions of each of the beams. Figure 2(d) shows the probability distribution of the centre of mass of the cell for a total simulation time of $5\rm {s}$, while Fig. 2(e) shows the orientation of the cell with respect to the fixed reference frame as a function of the simulation time. It is important mentioning that in the current configuration a rotation around the y-axis ($\beta$) would be a rotation around the cell axis of symmetry and therefore completely irrelevant. By extracting the average values for each degree of freedom, it is possible to compare the final equilibrium configuration obtained with the NN and with GO. Indeed, the average values obtained with the predictions of the NN and GO methods agree well and are also in good agreement with previously reported values [21], Table 1.

Table 1. Equilibrium position and orientation for a RBC in a double-beam OT as found with geometrical optics (GO) and with neural networks (NN). For each parameter we report the average and the standard deviation.

View Table

Moreover, the biggest advantage of using the NN for numerical simulations is a consistent decrease in the simulation time required to achieve the same precision (the NN is two orders of magnitude faster). Since the NN shows a higher computational efficiency, hereafter, we make use of the NN prediction to simulate the Brownian dynamics of an optically trapped RBC.

We therefore move to extract quantitative information on the trap constants. Initially we analyse the hydrodynamics of the RBC, since non-spherical particle could have an intrinsic roto-translation coupling due to their peculiar shape [37]. In our case, the diffusion tensor D does not show any strong roto-translation coupling; therefore, we do not expect to find any strong correlation in the cell’s motion intrinsically due to the RBC’s hydrodynamics. Still, optically trapped non-spherical particles could show roto-translation coupling in their motion as previously observed by others. In this framework, the normalised auto-correlation function (ACF) has been successfully used to extract quantitative information about the trapping constants [41,42].

We first evaluate the spatial ACFs $\left (C_{\rm {xx}}\left (\tau \right ), C_{\rm {yy}}\left (\tau \right ), C_{\rm {zz}}\left (\tau \right )\right )$ of the particle centre of mass trajectories. $C_{\rm {xx}}\left (\tau \right )$ and $C_{\rm {zz}}\left (\tau \right )$ decay as a single exponential with characteristic decay frequencies $\omega _{\rm {x}}=28{\ \mathrm {s}}^{-1}$ and $\omega _{\rm {z}}=6.4\ \mathrm {s}^{-1}$. Contrariwise, $C_{\rm {yy}}\left (\tau \right )$ is well fitted with a double exponential with characteristic frequencies $\omega _{y,1}=42{\ \mathrm {s}}^{-1}$ and $\omega _{y,2}=2.7\ \mathrm {s}^{-1}$, Fig. 3(a). We associate the fast decay rate to the translation, while the slower decay can be related to rotation around the $x$-axis ($\alpha$) induced by a motion along the y-direction. The values of the normalized cross-correlation function between $\alpha$ and y at zero time lag ($C_{\alpha y} (0) = -0.368$) further confirm a roto-translation coupling, Fig. 3(c). [43] Fig. 3(b) shows a density plot of the rotation around the $x$-axis ($\alpha$) as function of the motion along the $y$-direction. Here it can be seen a moderate negative correlation which suggests that the RBC rotates as it moves away from $y_{\rm {eq,2}}$, and undergoes to an “oscillating” motion about the equilibrium configuration where it is stably confined. To better comprehend this correlation we simulate $F_{\rm {y}}\left (\alpha \right )$ (Fig. 3(d)) and $\tau _{\rm {x}}\left (y\right )$ (Fig. 3(e)) which undoubtedly shows the coupling between the motion along $y$ and $\alpha$. Actually, the cell in its “folded” position (i.e. $\alpha = 90^\circ$) is constantly subjected to a force along the $y$-direction that moves the particle away from $y_{\rm {eq,2}}$ which in turns induces a rotation around the $x$-direction. On the other hand, as extensively described by Tognato et al., the transverse forces and torques components confine the cell in its “folded” configuration [21]. The overall consequence of these stable and unstable equilibria is a “circulating motion” of the cell within the optical trap about the equilibrium configuration. This would suggest that the coupling is intrinsically due to particle shape and to the optical trap rather than to the hydrodynamic of the particle.

 

Fig. 3. a) Translational autocorrelation function. The solid lines are exponential fits. $C_{xx}\left (t\right )$, $C_{zz}\left (t\right )$, decay as single exponential while $C_{yy}\left (t\right )$ as double exponential. b) $y-\alpha$ correlation shown as density plot. c) Normalised cross-correlation function between the rotation around the x-axis ($\alpha$) and the y-displacement (red line exponential fit). Both $F_y\left (\alpha \right )$ d) and $T_x\left (y\right )$ e) reveal unstable equilibrium when the cell is tilted of 90$^{\circ }$ around the x-axis (i.e. RBC in its folded position)

Download Full Size | PDF

Nevertheless we do not consider directly the membrane properties, any pathology that affects the morphology of a RBC can be monitored by our method. In fact, Brownian motion can be particularly useful in detecting shape asymmetry [44]. Shape irregularities in the RBC morphology could be caused by a pathogen infection as well as genetic inherited diseases. These shape deformities are mainly due to a drastic change in the biomechanical properties of RBC. For example, patients affected by diseases like malaria [10] or sickle cell disease [45] present erythrocytes with various level of morphological distortion. These anomalies could lead to a divergence from an ideal Brownian motion of a optically trapped healthy RBC which can be modeled with our methodology. For example, these roto-translation deviations from the ideal Brownian motion could be used in detecting pathogenic or genetic inherited diseases, helping in an early diagnose of disease.

Lastly, we extract average values and the standard deviations for the force constants ($k_{2,x}=\frac {\omega _{\rm {x}}k_{b}T}{D_{\rm {xx}}}=0.166\pm 0.024 \,\frac {\rm {pN}}{\mu \rm {m} \cdot \rm {mW}}$, $k_{2,y}=\frac {\omega _{\rm {y}}k_{b}T}{D_{\rm {yy}}}=0.218\pm 0.025 \,\frac {\rm {pN}}{\mu \rm {m} \cdot \rm {mW}}$, $k_{2,z}=\frac {\omega _{\rm {z}}k_{b}T}{D_{\rm {zz}}}=0.005\pm 0.001 \,\frac {\rm {pN}}{\mu \rm {m} \cdot \rm {mW}}$). These values are in excellent agreement with a previously reported work [21]. Similarly to the translational motion, we calculate $C_{\alpha \alpha }\left (\tau \right )$ and $C_{\gamma \gamma }\left (\tau \right )$. $C_{\alpha \alpha }\left (\tau \right )$ and $C_{\gamma \gamma }\left (\tau \right )$ decay as a single exponential and the respective trap constant are: $k_{\alpha }=\frac {\omega _{\alpha }k_{b}T}{D_{\rm {\alpha \alpha }}}=0.352\pm 0.096 \,\frac {\rm {pN}\cdot \mu \rm {m}}{\rm {rad}\cdot \rm {mW}}$ and $k_{\gamma }=\frac {\omega _{\gamma }k_{b}T}{D_{\rm {\gamma \gamma }}}=1.587\pm 0.382 \,\frac {\rm {pN}\cdot \mu \rm {m}}{\rm {rad}\cdot \rm {mW}}$. We do not analyse the dynamics around $\beta$ since the cell is not confined about this axis.

3.3 Triple beam optical tweezers

As previously suggested, one of the greatest advantage of using a NN instead of GO is the significant lowered computation time, especially when a very high number of light rays is needed (e.g. a triple- or four beams optical tweezer). Now, we exploit this feature to investigate the equilibrium orientation and position of a RBC with a reconfigurable triple-beam OT.

If directly trapped, a healthy biconcave RBC can assume two different and alternative orientations within the optical trap depending on the number of beams used for trapping [14,21]. In a double-beam OT, the major axis of a RBC is parallel to the optical axis and the beam foci are contained in the cell, known as “folded” configuration [11]. On the contrary, if three or four beams arranged in symmetric configurations are used (i.e. beams foci on the vertex of equilateral triangle or a square), the major axis of the cell is confined to be orthogonal to the optical axis (i.e. $\alpha = 0 ^\circ$), configuration referred to as ’flat’ configuration [46]. Here, we sought for alternative (and intermediate) RBC equilibrium configurations in respect to the well-known “folded” and “flat” ones.

We consider a trap configuration that is intermediate to those traps able to trap the cell in its “folded” or “flat” configuration. We consider a triple-beam optical tweezers (TBOT) composed by three identical and tightly focused Gaussian laser beams. Two beams are always arranged along the $x$-axis in a diametrically opposite location on the thickest portion of the cell (white crosses in Fig. 4(a)). A third beam (yellow cross in Fig. 4(a)) can be translated over the thickest portion of the cell and is used to counteract $T_{\rm {x}}$ generated by the two fixed beams. For simplicity, henceforth, the position of the moving beam is described by a polar co-ordinates system in the $x-y$ plane. Its location is defined by a single angle ($\zeta$), and the distance from the origin is fixed and equal to the radius of the thickest portion of the cell ($2.76 \mu \rm {m}$), Fig. 4(a).

 

Fig. 4. a) Schematic depiction of the triple-beam optical trap and the polar co-ordinates system used to identify the position of the moving beam, and respective beam definition as beam 1, 2 and 3. (b) Force-field acting on the RBC located on a grid of coordinates in the x-y plane for $\zeta = 45^\circ$. The color scale indicates the magnitude of the total force acting on the x-y plane, while the grey arrows indicate the direction of the force. (c) Three-dimensional trajectories of the cell centre of mass over a simulation time of 10 s for different $\zeta$, and the average values for the last second of simulation. d) Polar ($\phi$) and e) azimuthal ($\theta$) orientation of the RBC as a function of the simulation time. Average orientations are measured over the last second of the simulation. The error bar represents the standard deviation. f,g) Final equilibrium configuration for a RBC in the reconfigurable triple beam optical trap for $\zeta =15^\circ$ and $90^\circ$ respectively. The blue dot indicates the center of mass of the RBC while the red stars indicates the beams’ foci. The numbers indicate the position of each of the beams.

Download Full Size | PDF

Next, we proceed with the identification of the positional and translational equilibria. As a first step in our investigation, we simulate a force-field acting on the cell for $\zeta = 45^\circ$ to appreciate the effect of the potential landscape on the RBC. In this simulation, the cell is in its “flat” configuration and located at $z = 0$. It can be seen that the light pattern creates a very complex force-field (Fig. 4(b)). Non-negligible optical forces act simultaneously along the $x-$ and $y$-direction for every location of the cell. The complexity of the force-field makes it extremely difficult to identify the equilibrium positions (i.e., point in space where a specific force component vanishes with negative slope). This process would require several reiterations for every degree of freedom, rendering the process labour intensive. However, noting that if a particle is subjected to an optical potential it falls into the equilibrium position/orientation, it would be possible to identify the equilibrium configuration studying its dynamics as suggested by Cao et al. [47]. From symmetry arguments, the effect of different locations of beam 3 can be understood as restricting $\zeta$ in the interval $[0^\circ,90^\circ ]$ as schematically depicted in Fig. 4(a). Moreover, since we are looking for alternative equilibrium configurations (or to a transition from a “flat-like” to “folded-like” configuration), it is also rational to disregard every position where two beams are too close to each other (i.e. $\zeta < 15^\circ$), which should induce a “folded” configuration. Thus, the position of beam 3 can be restricted to $15^\circ \leq \zeta \leq 90^\circ$. To evaluate the effect of the reconfigurable optical trap, $\zeta$ is sampled every $15^\circ$, and for each $\zeta$ the Brownian dynamics are simulated for a 10 s trajectory starting from a RBC positioned in its ’flat’ configuration ($\theta = 0^\circ$ and $\phi = 0^\circ$) centred at $(0,0,0)$. The simulation finishes once the cell equilibrates around a stable position and orientation. The final position and orientation are then given as the average position and orientation with the standard deviation of the last second of the simulation. Figure 4(c) shows the 3D trajectories of the RBC’s CM obtained from the simulations carried out for different $\zeta$. Here, while $x$ and $y$ equilibrium positions remain close to the origin for different angles, the equilibrium in $z$ does depend on $\zeta$. In particular, for $\zeta < 30^\circ$, $z_{eq}< 0 \mu \rm {m}$ and for $\zeta > 45^\circ$, $z_{eq}> 0 \mu \rm {m}$, Fig. 4(c). We anticipate that for $\zeta \leq 30^\circ$, the cell is in its “folded” configuration, Fig. 4(d) and Fig. 4(f). This is due to a combination of the light intensity distribution and the cell configuration within the trap. In fact, when the cell is in its “folded” position, the cell’s major axes are parallel to the direction of propagation of the light beam. In this condition, more highly converging “light rays” strike the biggest faces of the RBC. This increases significantly the gradient force ($F_g$). Simultaneously, while in folded position, the scattering force ($F_s$) decreases appreciably because of the smaller geometrical cross-section of the cell. However, if $\zeta$ increases, this effect is less pronounced since the light rays strike the cell less symmetrically, and for $\zeta = 30^\circ$, $z_{eq} \sim -0.2 \mu \rm {m}$. Conversely, for $\zeta \geq 45^\circ$ a net shift in the axial position is evidenced ($z_{eq} \sim 0.8 \mu \rm {m}$), and this is due to a sequential shifting from the “folded-like” configuration to a “flat-like” configuration, Fig. 4(c) and Fig. 4(d). Much more interesting is the analysis of the rotational equilibrium. In Fig. 4(d) are shown the polar orientation ($\phi$) of the cell as a function of the simulation time for different locations of the moving beam (i.e. various $\zeta$). It is evident that $\zeta$ strongly influences the final polar orientation of the cell, Fig. 4(d). In particular, as beam 3 approaches beam 2, the cell tilts more until it reaches the “folded” configuration (i.e. $\phi = 90^\circ$) for $\zeta = 30^\circ$. Analysing the final orientation of the cell in more detail, it is possible to discriminate between three different regions. When the two beams are close to each other, $\zeta \leq 30^\circ$, the cell is in the “folded” configuration. If $30^\circ \leq \zeta \leq 75^\circ$, the RBC’s tilting seems to vary linearly with $\zeta$, from a “folded-like” configuration to “flat-like” configuration. The last region is for $\zeta \geq 75 ^\circ$, where the cell tilting cannot be decreased further, Fig. 4(d). It is also interesting to note the minor effect that $\zeta$ has on $\theta$. Here, we define $\theta = 0^\circ$ when $z_{p}$ (defined in Fig. 1(a)) is pointing along the positive x-direction. For example, in the simple case of a double beam optical tweezers, the cell plane point towards the positive ($\theta = 90^\circ$) or negative ($\theta = 90^\circ$) y-axis. Either direction are equally plausible due to the symmetry of the cell. Therefore, in the case of a triple beam optical tweezers, the induced cell rotation around the z-axis is relatively small for different location of beam 3. In fact, the cell rotates at most of $\approx 10^\circ$. Yet, the rotation can be explained with the tendency of maximizing the overlapping volume between the trapped particle and the illuminating beam in order to minimize the energy of the system. This can be well understood in the discrete dipole approximation. In particular, for $\zeta =45 ^\circ$ it is possible to obtain the highest cell’s tilting around the z-axis, Fig. 4(e). For every other $\zeta$, the tilting of the RBC around the z-direction decreases towards $\theta =90^\circ$.

4. Conclusions

Although geometrical optics enables the calculation of optical forces exerted on a trapped RBC, achieving an accurate calculation had traditionally been a task known for its time-consuming nature. In this work we have demonstrated that by using NN, one can significantly increase the speed of calculation without compromising the accuracy. This enhancement of the force calculations has allowed us to explore systems that were almost impossible to tackle with the conventional method for optical force calculation. In particular, we have focused on the analysis of the dynamics of trapped RBC with multiple beams. This can potentially allow determination of the best trapping configuration and to minimize the incident laser power and therefore reduce the risk of photodamaging the trapped cells. Importantly, the proposed strategy can be readily extended to investigate position and orientational control over other complex-shaped particles using different beam configurations as well, broadening its applicability and potential impact.