CN109190558B - Method for monitoring real-time three-dimensional dynamic behavior of particles - Google Patents

Abstract

The invention discloses a monitoring method for the real-time three-dimensional dynamic behavior of particles. A coaxial digital holographic microscope is used to record the real-time hologram, the scattered light field of the particles is numerically reconstructed, and the local maximum light intensity is used to obtain the three-dimensional position of the particles; The three-dimensional positions are connected to form a trajectory; the root mean square end-to-end distance of the trajectory under different time intervals is calculated to obtain the MSD-Δt curve; the curve is fitted to classify the trajectory; the instantaneous velocity of each particle after classification, the instantaneous velocity vector and the z-up direction are calculated. The included angle and the average velocity of each trajectory; determine the position of the interface, calculate the axial distance from each trajectory point to the interface, the spatial distribution of particle density and velocity, and the distribution of θ; the invention simultaneously measures multiple particles, and the axial positioning accuracy It can obtain real-time three-dimensional coordinates of particles; it does not require fluorescent labels, and it is non-destructive imaging; the imaging depth reaches hundreds of microns, which is very suitable for dynamic three-dimensional tracking and analysis of the activities of various particles.

Description

A method for monitoring real-time three-dimensional dynamic behavior of particles

technical field

The invention relates to the research field of three-dimensional position monitoring of particles, in particular to a method for monitoring real-time three-dimensional dynamic behavior of particles.

Background technique

Microparticles refer to particles ranging in size from submicron to hundreds of microns, including bacteria, fungi, viruses, nanoparticles, colloids and cells. Real-time monitoring of the three-dimensional dynamic behavior of particles in the medium and near the interface can be used to study the response of the particles to the environment (temperature, electric field, pH, etc.) and the interaction between the particles and the interface, and reveal the mechanism of such processes.

Traditional optical microscopes can only record the two-dimensional coordinates of the particles at the focal plane, and are not suitable for tracking their three-dimensional motion. The confocal microscope can record the information of particles at different heights by means of layer scanning, and realize three-dimensional tracking. However, the time resolution of confocal microscopy is greatly limited by the scanning speed of the instrument, and it is not suitable for real-time tracking of mobile and flexible particles such as bacteria. Total internal reflection microscopy utilizes the exponential decay of evanescent wave intensity in the axial direction, which can monitor the three-dimensional movement of particles in real time and non-destructively, and the axial positioning accuracy is sub-nanometer to several nanometers. However, its imaging range in the axial direction is only a few hundred nanometers, and the spatial variation range of particle activity can reach tens to hundreds of micrometers, which limits the application of this technology in three-dimensional tracking of particles.

SUMMARY OF THE INVENTION

The purpose of the present invention is to overcome the shortcomings and deficiencies of the prior art, and to provide a method for monitoring the real-time three-dimensional dynamic behavior of particles.

The object of the present invention is achieved through the following technical solutions:

A method for monitoring real-time three-dimensional dynamic behavior of particles, comprising the following steps:

S1. Record the holograms at different times on the sample to obtain the original holographic image, record the background image in the blank area of the sample, or average the light intensity of the original holographic image to obtain the background image, and the sample is a particle;

S2. Perform background subtraction on the original holographic image through the background image to eliminate background noise;

S3, performing numerical reconstruction on the background-deducted holographic image to obtain the three-dimensional reconstructed intensity distribution information of each particle at different times;

S4. Reconstructing the three-dimensional intensity distribution information of the particles, filtering through a threshold value and finding the local maximum value of the light intensity, to obtain the candidate three-dimensional positions of each particle at different times;

S5. Set the displacement threshold, determine the three-dimensional positions of the same particle at different times, and connect them to form a three-dimensional trajectory; calculate the root mean square end distance MSD of the trajectory under different time intervals Δt, and obtain the MSD-Δt curve;

S6. Fit the MSD-Δt curve, and classify the trajectory of the particles according to the motion mode by the index obtained by fitting; calculate the instantaneous velocity V of each trajectory point and the clip between the instantaneous velocity vector and the upward direction of z in different modes. Angle θ, the average speed V _T of each track is the average value of the instantaneous speed of all track points in the track;

where the instantaneous speed is:

Among them, i is the frame number of the track point in the track, i=2,3...N, N is the total number of frames of the track; (x _i , y _i , z _i ) and (x _i-1 , y _{i- 1} , z _i-1 ) is the position of the same particle in two consecutive frames, dt is the time interval of two consecutive frames;

The angle θ between the instantaneous velocity vector and the z-up direction:

为z朝上方向的方向向量，

为速度矢量，

|a|为

的模，|b|为

的模；in,

is the direction vector in the upward direction of z,

is the velocity vector,

|a| is

the modulus of , |b| is

the model;

S7. Take the lowest axial position among all the trajectory points as the position z ₀ where the interface is located, and the axial position of each trajectory point is _zi ; calculate the axial distance z _s from each trajectory point to the interface, z _s = _zi -z ₀ ;

S8. Calculate the spatial distribution of the density and velocity of the particle and the distribution of θ under different motion modes; analyze the three-dimensional dynamic behavior of the particle and the interaction between the particle and the interface through the motion characteristic quantity.

Further, the average light intensity of the original holographic image is specifically calculated by the formula as:

Among them, I _b (x, y) is the gray value of the pixel at the position (x, y) in the background image, N is the total number of frames of the holographic image, t is the time, and I _t (x, y) is the original image at time t. The gray value of the pixel at the (x, y) position in the holographic image;

Further, in step S2, the background subtraction of the original holographic image is specifically: reading the gray value of each pixel in the original holographic image and the background image, and the gray value of each pixel in the holographic image after background subtraction. The value Is ( _x ,y) is:

I _s (x,y)=I _t (x,y)-I _b (x,y),

Among them, I _b (x, y) is the gray value of the pixel at the position (x, y) in the background image, and I _t (x, y) is the gray value of the pixel at the position (x, y) in the original hologram at time t. degree value;

Further, in the step S3, the specific process is:

U(r,z)=FT ^-1 ( _FT (Is(r,0)·H(q,-z))),

Among them, h(r,-z) is the propagation operator, r is the initial lateral coordinate of the particle, z is the initial axial coordinate of the particle; i is the imaginary unit; _k is the wave number; R is the light propagation distance; Is is the sample The light intensity of ; FT ^-1 is the inverse Fourier transform; FT is the Fourier transform; H(q,-z) is the Fourier transform of h(r,-z);

The value of the reconstructed axial interval is greater than the imaging range of the instrument in the axial direction, and the step is the size of the pixel divided by the magnification of the objective lens;

Further, the threshold filtering, the light intensity threshold is less than 10%;

Further, in the step S5, the specific process is: according to the movement characteristics of the particles, a displacement threshold value T is set, and if the displacement value between two positions in two consecutive frames is less than the threshold value, the two positions belong to the same particle; at the same time, The length of the set trajectory should not be less than W points, W=30, to ensure the correctness of the trajectory connection and the root mean square end distance calculated in the subsequent steps;

Further, the displacement threshold T is obtained by multiplying the particle velocity V _s by the exposure time t _e , namely: T=V _s t _e ;

Further, the fitting of the MSD-Δt curve is specifically:

The points in the first 10% of the MSD-Δt curve are fitted by the following formula:

MSD(Δt)=6D(Δt) ^ν ,

The index ν obtained by the fitting divides the trajectory of the particle according to the motion mode, and the result is: the index ν<0.95 is the first mode; the index 0.95≤ν≤1.05 is the second mode; the index 1.05<ν is the third mode;

Further, the first mode is a restricted diffusion mode; the second mode is a free diffusion mode; the third mode is an active motion mode;

Compared with the prior art, the present invention has the following advantages:

The invention uses the Rayleigh-Sommerfeld algorithm to locate the three-dimensional position of the particles, can measure multiple particles at the same time, and the axial positioning accuracy reaches the sub-micron level; the real-time three-dimensional coordinates of the particles can be obtained; it does not need fluorescent markers, and is a non-destructive imaging method . Its imaging depth can reach tens of hundreds of microns, which is very suitable for dynamic tracking of the activities of various types of particles.

Description of drawings

Fig. 1 is a method flow diagram of a method for monitoring the real-time three-dimensional dynamic behavior of particles according to the present invention;

2 is a structural block diagram of a coaxial digital holographic microscope imaging system in the embodiment of the present invention;

Figure 3 is the original hologram of Escherichia coli in the embodiment of the present invention;

FIG. 4 is a background image obtained by averaging the light intensity of a holographic image sequence in the embodiment of the present invention;

Fig. 5 is the hologram deducted background in the described embodiment of the present invention;

Fig. 6 is the movement track diagram of Escherichia coli near the interface in the embodiment of the present invention;

Fig. 7 is the MSD-Δt curve diagram of the restricted diffusion trajectory in the described embodiment of the present invention;

Fig. 8 is the MSD-Δt curve diagram of the active motion trajectory in the embodiment of the present invention;

Fig. 9 is the density distribution diagram of Escherichia coli in active movement mode in the embodiment of the present invention;

10 is a spatial distribution diagram of Escherichia coli in active movement mode in the embodiment of the present invention;

11 is a distribution diagram of the included angle θ between the velocity vector of Escherichia coli in the active movement mode and the upward z direction in the embodiment of the present invention.

Detailed ways

The present invention will be described in further detail below with reference to the embodiments and the accompanying drawings, but the embodiments of the present invention are not limited thereto.

Example

A method for monitoring real-time three-dimensional dynamic behavior of particles, as shown in Figure 1, includes the following steps:

The first step: use the coaxial digital holographic microscope imaging system to record holograms at different times on the sample to obtain the original holographic image, record the background image in the blank area of the sample, or average the light intensity of the original holographic image to obtain the background image; the The sample is a particle; the calculation process of the average light intensity is as follows:

Among them, I _b (x, y) is the gray value of the pixel at the position (x, y) in the background image, N is the total number of frames of the original holographic image, t is the time, and I _t (x, y) is the time t The gray value of the pixel at the (x, y) position in the original holographic image;

In this embodiment, the microparticles are selected as wild-type Escherichia coli;

The coaxial digital holographic microscope imaging system is shown in FIG. 2 , including two mutually perpendicular LED light sources, a first light source 1 and a second light source 1 ′, a first plane mirror 7 , a sample stage 8 , and a microscope objective lens 9 , a second plane reflection 10 and a sCMOS camera 11; wherein, the two LED light sources are designed with a coaxial cage light path, and the light path is simple and flexible to change the incident wavelength; the LED light sources both include a first LED light 2 and a second LED Lamp 2', first microscope objective lens spatial filter 3 and second microscope objective lens spatial filter 3', first diaphragm 4 and second diaphragm 4', and first convex lens 5 and second lens with focal length f=50mm Convex lens 5', LED light source 1 also includes a third plane mirror 6; microscope objective lens (10, 20, 40, 60 and 100 times magnification) is installed on the nosepiece, which can flexibly switch the objective lens.

Step 2: Perform background subtraction on the original holographic image through the background image to eliminate background noise; specifically: read the gray value of each pixel point in the original holographic image and the background image, and each pixel point in the holographic image after background subtraction The gray value of Is ( _x , y) is:

I _s (x,y)=I _t (x,y)-I _b (x,y),

Among them, I _b (x, y) is the gray value of the pixel at the position (x, y) in the background image, and I _t (x, y) is the gray value of the pixel at the position (x, y) in the original hologram at time t. degree value; FIG. 3 is the original hologram of Escherichia coli in the embodiment of the present invention; FIG. 4 is the background image obtained by averaging the light intensity of the holographic image sequence in the embodiment of the present invention; The background-subtracted hologram in the above embodiment.

Step 3: According to the Rayleigh-Sommerfeld algorithm, numerically reconstruct the holographic image deducted from the background, and obtain the intensity distribution information of the reconstructed light field of all E. coli within the field of view in the set three-dimensional space. The calculation is as follows:

U(r,z)=FT ^-1 (FT(Is(r,0)·H(q,-z))),

Among them, h(r,-z) is the propagation operator, r is the initial lateral coordinate of E. coli, z is the initial axial coordinate of E. coli; i is the imaginary unit; k is the wave number; R is the light propagation distance; I _s is the light intensity of the sample; FT ^-1 is the inverse Fourier transform; FT is the Fourier transform; H(q,-z) is the Fourier transform of h(r,-z);

Step 4: For the three-dimensional reconstruction intensity distribution information of Escherichia coli obtained above, obtain the candidate three-dimensional positions of each particle at different times through threshold filtering and finding the local maximum value of the light intensity; the specific steps are to set a light intensity first. The threshold, usually less than 10%, is set to 7% here to filter out the noise in the reconstructed light field, and then the edge length is w (taking bacteria as an example, w takes an integer value close to the sum of the bacterial flagella and body length) Find the local maximum light intensity point by point in the cube. Here, in the cube with a side length of 13 μm, find the three-dimensional position of the point of maximum light intensity point by point and record it.

Step 5: Set the displacement threshold to 3.7 μm, determine the three-dimensional position of the same E. coli at different times, and connect them into a three-dimensional trajectory; if the displacement value between the two positions in two consecutive frames is less than the threshold, then the two positions Belong to the same Escherichia coli; the track length is set to not less than 30 points to ensure the correctness of the root mean square end distance calculated by the track connection and subsequent steps; Fig. 6 is the movement of Escherichia coli near the interface in the embodiment of the present invention Trajectory diagram; calculate the RMSD distance from the end of the trajectory to MSD at different time intervals Δt, and obtain the MSD-Δt curve:

MSD(Δt)=<|r(t ₀ +Δt)-r(t ₀ )| ² >,

Among them, the value of Δt is 0.05ns, n=1,2,3...N-1, and N is the length of the trajectory; Fig. 7 is the MSD-Δt curve diagram of the restricted diffusion trajectory in the embodiment of the present invention; Fig. 8 is the MSD-Δt curve diagram of the active motion trajectory in the embodiment of the present invention;

Step 6: Fit the MSD-Δt curve and fit the first 10% points of the curve:

MSD(Δt)=6D(Δt) ^ν ,

The index obtained by fitting classifies the trajectory of Escherichia coli according to the movement mode, which is divided into restricted diffusion mode, free diffusion mode, and active movement mode; among them, the index ν＜0.95 is the restricted diffusion mode; the index 0.95≤ν≤1.05 Free diffusion mode; index 1.05<ν is active motion mode; this example aims to analyze the complete three-dimensional dynamic behavior of Escherichia coli near the interface, trajectories with ν≥1 are regarded as active motion, and subsequent data processing is limited to active motion. Escherichia coli;

为该轨迹中所有轨迹点瞬时速度的平均值；Calculate the instantaneous velocity V of each trajectory point and the angle θ between the instantaneous velocity vector and the upward direction of z in different modes, and the average velocity of each trajectory

is the average value of the instantaneous speed of all trajectory points in the trajectory;

where the instantaneous speed is:

Among them, i is the frame number of the track point in the track, i=2,3...N, N is the total number of frames of the track; (x _i , y _i , z _i ) and (x _i-1 , y _{i- 1} , z _i-1 ) is the position of the same particle in two consecutive frames, dt is the time interval of two consecutive frames;

The angle θ between the instantaneous velocity vector and the z-up direction:

为z朝上方向的方向向量，

为速度矢量，

|a|为

的模，|b|为

的模in,

is the direction vector in the upward direction of z,

is the velocity vector,

|a| is

the modulus of , |b| is

the mold

Step 7: Take the lowest axial position of all the trajectory points as the interface position z ₀ , and the axial position of each trajectory point as _zi ; calculate the axial distance z _s from each trajectory point to the interface, namely: z _s = z _i -z ₀ ;

Step 8: Calculate the spatial distribution of the density and velocity of the particles and the distribution of θ under different motion modes; analyze the three-dimensional dynamic behavior of the particles and the interaction between the particles and the interface through the motion characteristics;

Here, take dz=5μm (much greater than the axial positioning resolution of 400nm), dθ=5°, and calculate the spatial distribution of the density and velocity of the actively moving Escherichia coli and the distribution of θ;

The three-dimensional dynamic behavior of Escherichia coli and the interaction between Escherichia coli and the interface were analyzed by the above motion characteristic quantities. It can be seen from the figure that the movement of Escherichia coli near the interface is mainly affected by hydrodynamics. The specific analysis is as follows: The spatial distribution of the density of E. coli near the interface shown in Figure 9 conforms to the theoretical model of particle distribution near the interface under the action of hydrodynamics:

In the formula, L _⊥ is the action distance of the hydrodynamic action, H is the thickness of the sample cell, p is the force dipole, η is the viscosity of the medium, and D is the diffusion coefficient; as shown in Figure 10, when z<5μm, E. coli moves The phenomenon that the velocity increases with the decrease of z and the phenomenon that θ is mainly distributed at 90° in Fig. 11 is consistent with the theoretical model that the particles are subjected to hydrodynamic action near the interface and tend to move in the direction parallel to the interface, and the velocity increases.

The above-mentioned embodiments are preferred embodiments of the present invention, but the embodiments of the present invention are not limited by the above-mentioned embodiments, and any other changes, modifications, substitutions, combinations, The simplification should be equivalent replacement manners, which are all included in the protection scope of the present invention.

Claims (9)

1. A method for monitoring real-time three-dimensional dynamic behavior of particles is characterized by comprising the following steps:

s1, recording holograms of different moments on a sample to obtain an original holographic image, and recording a background image in a blank area of the sample, or carrying out light intensity averaging on the original holographic image to obtain a background image, wherein the sample is a particle;

s2, carrying out background subtraction on the original holographic image through the background image to eliminate background noise;

s3, carrying out numerical reconstruction on the background-subtracted holographic image to obtain three-dimensional reconstruction intensity distribution information of each particle at different moments;

s4, carrying out three-dimensional reconstruction on intensity distribution information of the particles, and filtering and searching a local maximum value of light intensity through a light intensity threshold value to obtain alternative three-dimensional positions of the particles at different moments;

s5, setting a displacement threshold, determining the three-dimensional positions of the same particle at different moments, and connecting the three-dimensional positions to form a three-dimensional track; calculating the root-mean-square end distance MSD of the track at different time intervals delta t to obtain an MSD-delta t curve;

s6, fitting the MSD-delta t curve, and classifying the particle track according to the motion mode by the index obtained by fitting; calculating the instantaneous speed V and the angle theta between the instantaneous speed vector and the upward z direction of each track point in different modes, and calculating the angle theta between each track point and the upward z directionAverage velocity

The average value of the instantaneous speed of all track points in the track is obtained;

wherein the instantaneous speed is:

wherein i is the frame number of the track point in the track, i is 2,3 … N, and N is the total frame number of the track; (x)_i，y_i，z_i) And (x)_i-1，y_i-1，z_i-1) Dt is the position of the same particle in two continuous frames, and dt is the time interval of the two continuous frames;

angle θ of instantaneous velocity vector to z-up direction:

wherein,

is a direction vector of the z-up direction,

in the form of a velocity vector, the velocity vector,

a is

Modulo b is

The mold of (4);

s7, taking the lowest axial position in all track points as the position z of the interface₀The axial position of each track point is z_i(ii) a Calculating the axial distance z between each track point and the interface_s，z_s＝z_i-z₀；

S8, calculating the spatial distribution of the density and the speed of the particles and the distribution of theta in different motion modes; and analyzing the three-dimensional dynamic behavior of the particles and the interaction between the particles and the interface through the motion characteristic quantity.

2. The method for monitoring the real-time three-dimensional dynamic behavior of particles as claimed in claim 1, wherein in step S1, the light intensity average of the original holographic image is calculated by the following formula:

wherein, I_b(x, y) is the gray value of the pixel at the (x, y) position in the background image, N is the total frame number of the original holographic image, t is the time, I_t(x, y) is the gray value of the pixel at the (x, y) position in the original holographic image at time t.

3. The method for monitoring the real-time three-dimensional dynamic behavior of particles as claimed in claim 1, wherein in step S2, the background subtraction is performed on the original holographic image, specifically: reading gray values of pixel points in the original holographic image and the background image, and subtracting the gray value I of the pixel points in the original holographic image after the background is deducted_s(x, y) is:

I_s(x,y)＝I_t(x,y)-I_b(x,y)，

wherein, I_b(x, y) is the gray scale value of the pixel at the (x, y) position in the background image, I_t(x, y) is the gray value of the pixel at the (x, y) position in the original hologram at time t.

4. The method for monitoring the real-time three-dimensional dynamic behavior of particles as claimed in claim 1, wherein the step S3 comprises the following steps:

U(r,z)＝FT^-1(FT(I_s(r,0)·H(q,-z)))，

wherein h (r, -z) is a propagation operator, r is an initial transverse coordinate of the particle, and z is an initial axial coordinate of the particle; i is an imaginary unit; k is the wave number; r is the light propagation distance; i is_sIs the light intensity of the sample; FT^-1Is inverse Fourier transform; FT is Fourier transform; h (q, -z) is the Fourier transform of H (r, -z);

the value of the reconstruction axial interval is larger than the axial imaging range of the instrument, and the stepping is that the size of the pixel is divided by the magnification of the objective lens.

5. The method as claimed in claim 1, wherein in step S4, the light intensity threshold is less than 10%.

6. The method for monitoring the real-time three-dimensional dynamic behavior of particles as claimed in claim 1, wherein the step S5 comprises the following steps: setting a displacement threshold T according to the motion characteristics of the particles, wherein if the displacement value between two positions in two continuous frames is smaller than the threshold, the two positions belong to the same particle; meanwhile, the length of the set track is not less than W points, and W is 30.

7. The method as claimed in claim 6, wherein the displacement threshold T is a particle velocity V_sMultiplied by the exposure time t_eThe result is: t ═ V_st_e。

8. The method for monitoring the real-time three-dimensional dynamic behavior of particles according to claim 1, wherein in step S6, the fitting of the MSD- Δ t curve specifically comprises:

the top 10% point of the MSD- Δ t curve was fitted by the formula:

MSD(Δt)＝6D(Δt)^ν，

dividing the track of the particles according to the motion mode by the index nu obtained by fitting, wherein the index nu is less than 0.95 and is a first mode; the index nu is more than or equal to 0.95 and less than or equal to 1.05 and is a second mode; the index nu of 1.05 is the third mode.

9. The method for monitoring the real-time three-dimensional dynamic behavior of particles as claimed in claim 8, wherein the first mode is a limited diffusion mode; the second mode is a free diffusion mode; the third mode is an active motion mode.

Images