A NEW ALGORITHM TO ANALYZE THE VIDEO DATA OF CELL CONTRACTIONS IN MICROFLUIDIC PLATFORMS

One of the rapidly developing trends in science is tissue engineering. Its goal is the creation of biological tissue with prescribed properties. One possible approach consists in culturing cells in microfluidic platforms (MPs) [1–4]. Such technologies allow for studying a proto-tissue or organ at the miscroscale resembling the principal features of the organ of. In general, MP can be successfully used to: – improve the technology of tissue cultivation; – determine the optimal conditions for creating tissues with given properties; – study the biophysical properties of interrelated groups of cells; – model pathological processes at the tissue level; – study the effect of external physical factors and therapeutic drugs on tissue; – develop and optimize clinical diagnostic tests by determining the individual susceptibility of the patient’s tissues to physical factors and therapeutic drugs. Various chemical biomarkers and mediums have been widely used to study the morphological and functional properties of tissues. For example, in [5, 6] medias and preparations for determining the proliferate activity of cells and phases of the cell’s life cycle, the number of dead and living cells, wave processes of spread of the front to stimulation with use and action potentials are given. The fundamental importance of using biomarkers is their specificity. However, we can identify a number of shortco mings of this approach. In most cases, biomarkers are quite expensive; their use requires special knowledge, skills, and appropriate equipment. A NEW ALGORITHM TO ANALYZE THE VIDEO DATA OF CELL CONTRACTIONS IN MICROFLUIDIC PLATFORMS


Introduction
One of the rapidly developing trends in science is tissue engineering.Its goal is the creation of biological tissue with prescribed properties.One possible approach consists in culturing cells in microfluidic platforms (MPs) [1][2][3][4].Such technologies allow for studying a proto-tissue or organ at the miscroscale resembling the principal features of the organ of.In general, MP can be successfully used to: improve the technology of tissue cultivation; -determine the optimal conditions for creating tissues with given properties; study the biophysical properties of interrelated groups of cells; model pathological processes at the tissue level; study the effect of external physical factors and therapeutic drugs on tissue; develop and optimize clinical diagnostic tests by determining the individual susceptibility of the patient's tissues to physical factors and therapeutic drugs.
Various chemical biomarkers and mediums have been widely used to study the morphological and functional properties of tissues.For example, in [5,6] medias and preparations for determining the proliferate activity of cells and phases of the cell's life cycle, the number of dead and living cells, wave processes of spread of the front to stimulation with use and action potentials are given.The fundamental importance of using biomarkers is their specificity.However, we can identify a number of shortco mings of this approach.In most cases, biomarkers are quite expensive; their use requires special knowledge, skills, and appropriate equipment.
When culturing cardiomyocytes on an MP, the most important functional test consists in the evalua tion of their contractility.One of the methods for evaluating mechanical contractions of cells is optical microscopy [7,2].The method does not require the use of additional biochemical preparations and the use of confocal microscopy.Cell contraction are evaluated on the basis of video data, which are obtained with conventional light microscopes.The technique allows for visualizing both spontaneous beats and evaluating cell oscillations under the influence of an external electric field, e.g., "pseudo-ECG" [2].To obtain quantitative characteristics of mechanical cell contractions without the use of auxi liary reagents, the task of computer processing and analysis of video images is posed.
In similar studies [2,8,9,10], where the amplitude of cardiomyocyte oscillations is estimated from the analysis of video images, the average change in brightness in the region of interest is measured.The size of this region usually is substantially smaller than the "useful area" of the microchip with cells and the method of estimating the oscillations is simple enough but has two significant drawbacks.First, the amplitude units of measurement is not calibrated, i.e., they are measured by levels of brightness.It means that it is impossible to estimate correctly cell oscillations in different areas of the microchip since the images are not always ideally illuminated by a uniform external light source, but also because of the variation in cell density in space.Analyzing the signal by luminance gradations (gray levels) when comparing the data of different experiments would require careful calibration.Another disadvantage is the distortion of the waveform characterizing the mechanical contraction of cells.This leads to a change in the balance of brightness (similar to the balance of mass or energy in physical tasks).In addition, the brightness of the optical image changes not only because of cells mechanical oscillations but also due to the polarization of the cell elements under the action of the external field.The overall brightness of the optical image of the cells varies nonlinearly with the change in the strength of the external field and the frequency of the effect.Thus, the known methods, for example, [2,8,9,10], register the very fact of cell contraction but substantially distort the shape and amplitude of the useful signal.
The aim of this study is the development of an algorithm for the analysis of the video data of mechanical contraction of cardiomyocytes on an MP to determine their functional and structural properties at the tissue level.

Materials and Methods
The developed algorithm for the analysis of video images was implemented in Matlab 2016.We analyzed the data of the contraction of cardiomyocytes cultured in an MP.The video data of cardiomyocyte contraction was recorded in the file format "mp4", where the color information in a pixel was recorded in the "RGB24" format.The recording rate was 23.83 frame/s.Subsequently, the images were transferred from the RGB format to grayscale images.
The geometric and functional features of the microchip described in [8].Cells are embedded in fibrin hydrogel and then injected into the microchips, cultured in DMEM high glucose [11] and stimulated according to the following 3 conditions: -Condition A: without any stimulation for 5 days.
-Condition B: without any stimulation for 3 days and subsequently electrically stimulated for the resting 2 days (day 4 and day 5).The electric pulse is 0,6 V amplitude (electric field 5 V/cm, current density 14.8 mA/cm 2 ), 2 ms positive, 2 ms negative and 0 V for the resting 996 ms; frequency -1 Hz.
-Condition C: without any stimulation for 3 days and subsequently stimulated for the resting 2 days (day 4 and day 5).The electric pulse is 3 V amplitude (electric field 25 V/cm, current density 74 mA/cm 2 ), 2 ms positive, 2 ms negative and 0 V for the resting 996 ms; frequency 1 Hz.
Video image were obtained after 5 days of culture.To stimulate the contraction, the chips were stimulated by a unipolar rectangular signal with a duty ratio of 2. The amplitude of the pulses was va ried in the range 0-7 V, the frequency was 0-5.5 Hz.
The algorithm for the analysis of cardiomyocyte contraction on an MP consists of the following sequence of actions: 1. Video data are loaded.2. Images (definition of the geometrical sizes of pixels) are calibrated.
3. A rectangular region of interest is defined.4. Binary images of local maxima in each frame of the video sequence are created.
5. Binary images of local maxima and binarization of the resulting image at the α-level of 0.5 is performed.
6.In the neighborhood of the obtained points in the spatial window of 9 pixels on each frame, the radius (R) and angle (θ) formed by the vector between the center of mass of the pixel brightness are calculated.
7. The changing rate of R and θ in time are calculated : ∆R/∆t and ∆θ/∆t.
8. The Region of Interest (ROI) along X and Y axis are stratified.9.The calculated data R and θ, ∆R/∆t and ∆θ/∆t in the strata for each time frame of the video sequence are averaged in space.
10.The obtained data, i.e., the shape of the oscillations in time, spectral Fourier analysis, correlation of oscillations between strata are analyzed.
The determination of the geometric dimensions of the pixels was carried out with respect to the 300 μm distance between the geometric centers of the columns shown in Fig. 1.In our case one pixel was equal 1.43 µm.
After ROI building, local maxima (BW max (t)) are searched for each video frame.It is assumed that these points determine the coordinates of cardiomyocytes, since the elements of the cytoskeleton have the largest optical density, which is visualized as the brightest (almost white) areas on the images.It should be kept in mind that during contraction cells change their geometric dimensions and move for a certain distance in space.In addition, a noise exists in the images.Therefore, to find the average equilibrium coordinate of the cells, all values of BW max (t) were summed up: where n t is a number of frames.The resulting B Σ image was normalized to the maximum value and binarized at 0.5 (B Σ0.5 ).The obtained points with the value of 1 were identified with the geometric centers of the cells.This technique allows avoiding the use of the alignment procedure for the non-uniform illuminating images and significantly reduces the computer time for analyzing the video stream, i.e., there is no need to process every pixel of images.We note that each cell may have several non-zero pixels B Σ0.5 .However, it should be emphasized that for the video images with a higher degree of spatial resolution, the thresholding technique may be different, for example, a relatively narrow range of brightness can be used.
The oscillation analysis of an individual cardiomyo cyte was performed in the neighborhood of nonzero pixels B Σ0.5 , which we conditionally call "active image points" (AIP).The neighborhood of the AIP was defined by a rectangular area of ±4 pi xels, that is, the window size was 9 pixels (6.3 µm).The size of the window was determined by two factors.First, given the resolution of video images, the distance between non-zero pixels B Σ0.5 was 5-20 pi xels.The second factor was the requirement that if the window size in the microchip area was increased, where the cells were a priori absent, the measured signal (for example, the average brightness) would change minimally.The maximum window size should not exceed the maximum cell size of cardiomyocytes, which is about 24 µm.The relative error of the amplitudes of the oscillations (va riable component) bet ween the data obtained with the windows 9, 11, 13, 15 did not exceed 5-10 % in our studies.In this regard, to minimize the machine resources, a window of 9 pixels was selected (see Fig. 2).
For the analysis of cell contraction, two main parameters were chosen: the radius and the angle formed by the vector between the geometric center of the window (the neighborhood of the non-zero point in image B Σ0.5 ) and the center of mass of the brightness of the window (Fig. 3).
When stimulating the cell oscillations with an external electrical signal, the electric field lines are directed along the X or Y axis.The spatial orientation of cardiomyocytes with respect to field lines can be considered as random one.In connection with this, it is also rational to estimate the components of the radius of oscillations R x and R y : The time derivatives ∆ t R, ∆ t Θ were estimated as follows: where ∆t was 1/23.83 s; n t is a frame number of the video sequence (time).Thus, the vector D = (i, formed, where i is a sequence number of the AIP, X i and Y i are the coordinates of the i-th point. To average the data, the ROI space was divided into strata (subranges) along the X and Y axis.The number of striations along the given axis N strX or N strY was determined by the H.A. Sturges formula where N X and N Y are the number of pixels in the ROI along the X and Y axes, respectively.The vectors D i , which coordinates X i , Y i correspond to the given stratum, are averaged.
According to the time variation of R i , Θ i it is possible to evaluate the contraction mode of the cells, and also distinguish low-amplitude oscillations from noise random components.
The subsequent analysis of the signals consisted of a correlation analysis of the oscillations of the mean values of R, R x , R y , Θ, ∆ t R, ∆ t Θ between strata, and also in the analysis of Fourier spectra.In the striations along the X axis, the components of the contraction R y were analyzed, and, conversely, in the strata along the Y axis, the average oscillations of R x were estimated.
The image processing method described above substantially reduces these drawbacks.In addition, the method makes it possible to estimate the magnitude of cell contraction in absolute units of μm, and the rate of oscillations is μm/s, and to decompose the mechanical oscillations of cells into components.
Shape of the pulses of mechanical oscillations of cardiomyocytes during their electrical stimulation are presented in Fig. 4. It can be seen that the form of the oscillations of the average brightness tends to rectangular.The oscillations of the radius R have a more complex form.Such a form of signal (Fig. 4b) is similar to the data in [6].
Among all AIP a sufficiently large number of points not belonging to the cells is present.This is largely due to the fact that the environment and micro-inclusions, including inactive and dead cells, also perform mechanical contractions.These points can significantly affect the results of averaging the data.It is possible to distinguish the AIPs belonging to living cells from the others when comparing the shape of the pulses of mechanical oscillations: AIP non-belonging cells are characterized by low-amplitude almost chaotic oscillations.However, comparing the forms of oscillatory processes between all AIP requires a lot of computer resources.Alternatively, the following method is proposed.For each AIP the obtained values of ∆ t R max we normalize with respect to their maximal value max (∆ t R maxi ), where i = 1,…, n (n is the total number of AIPs).Those points for which the normalized values ∆ t R max are less than a certain thre shold level β are excluded from the analysis.In Fig. 5 an example of an estimation of the change in the shape of mechanical oscillations of cells is shown as a function of the level β (0 < β < 1).
The data in Fig. 5 show that the amplitude of the oscillations increases with increasing β (to 0.5); the shape of the pulse for β > 0 is quite different from the data for β = 0 (when all AIP are taken into account).The use of β > 0.5 is not rational, since in these cases the amount of AIP sharply decreases (Fig. 6), and the average data become sensitive to single relatively high-amplitude fluctuations.

Results
To approximate the average data, the equation was chosen, where f is the stimulation frequency, a and λ are parameters.
An example of the approximation of the average values of oscillations as a function of the stimulation frequency under different voltage on electrodes is presented in Fig. 7.
The quantitative values of the approximation of the experimental data are presented in the Table , where R 2 is a square of the correlation coefficient between the approximation function and experimental data.
According to the Table, the following conclusions can be drowned out: with the increase of the stimulating frequency, the amplitude of mechanical contractions decreases monotonically; for cells stimulated with electrical stimuli (groups B and C), the parameters a and λ decrease (in comparison with the group A), which indicates a decrease in the excitability of cells and the ability to perform mechanical work.A similar conclusion can be drown out by comparing groups B and C: an increase in the strength of the external electric field during the cell culture entails a decrease in the functional capability of cells.In group C, at the stimulation voltage of 4.95 V, there is a certain deviation from the presented general conclusions: the parameter a exceeds the corresponding value at the voltage of 4.5 V.This observation can be explained by the fact that the measurements were carried out on different microchips and every biological system has its own individual characte ristics.The presented analysis technique allows estimating the wave processes of spread contractions as well (see Fig. 8).
The wave character (when present) is clearly noticeable on two-dimensional time histograms ∆ t R along the X and Y axes, see Fig. 9a.It can be seen that the component of the rate of cell contraction ∆ t R Y lags behind the right edge of the microchip compared to the origin (the left edge of the microchip).At the same time, cell contractions in strata along the Y axis occur simultaneously (Fig. 9b).
Active cardiomyocytes on a microchip are distributed non-uniform in space.Therefore, the shape of the impulses can be different in different strata.When the stimulation conditions change, the most active regions can shift.The Fig. 10 shows examples of different pulse shape changes R x and R y in one of the chips.The oscillations along R x and R y axis maybe depended on the way cells remodeled the matrix and re-organized themselves.
At the stimulation frequency of 1 Hz the maxi mum amplitude of oscillations was observed (see Fig. 10a).The change in the constant component is also noted.The shape of the pulses changes, when the frequency increases.This phenomenon is observed at frequencies higher than 2 Hz. Figure 10b demonstrates that sometimes there are oscillations with opposite phase of contraction.Probably the phase shift between R x and R y depends on the number of cellular contacts between the cells and their orientation in space.Oscillations are antiphase if there are few contacts and the cells do not depend on neighbors.The absence of a pronounced phase shift indicates that the cells are oriented with respect to the external field randomly and the direction of the oscillations depends on the neighbors.Nevertheless, it can be stated that along the direction of the lines of external field force the amplitude of the oscillation is usually more pronounced (red curves).
In addition, in most cases asymmetry of the amplitude of cell contraction near two different stimulating electrodes is observed.For example, Fig. 9b shows that the speed and the amplitude of oscillations are higher in the upper part than in the lower one.This is due to the redistribution of charges in the medium.From the methods of electro-physiotherapy (galvanization, electrophoresis, low-frequency methods of exposure) it is known that under the anode a more acidic environment is formed and the sensitivity of tissues decreases, Table : The approximation functions for the maximum rate of cell contraction (μm/s) cultivated under different conditions, depending on the frequency and the stimulation voltage of the oscillations

Group of microchips
Parameters a and λ of approximation functions under the cathode there is an increase in pH and an increase in the sensitivity of tissues [12].In turn, this indicates that when testing the electromechanical properties of cells in a microchip as a complete system, it is necessary or to increase the duty cycle of the exciting external pulses and/or use to a bipolar signal.This will make it possible to restore the spatial distribution of ions.Since the mobility of positive and negative ions is different, one should also limit the time of one measurement and make substantially larger pauses between measurements on a single chip.With an increase in the oscillation frequency of more than 2 Hz, a low-frequency modulation of the pulse shape is often observed (see Fig. 10d).
It should be noted that at this stage of the study the analysis of the time variation of the parameters Θ, ∆ t Θ did not provide additional useful information.These parameters correlated well with the amplitude characteristics but had a slightly larger noise component.However Θ and ∆ t Θ can be of interest for cluster space-time analysis.The spatial heterogeneity in the mechanical activity of cells can also be considered with the use of maps of the oscillation rate distribution (Fig. 11).The size of one pixel on these maps corresponds to the size of the strata along the corresponding axis.As a rule, with increasing the stimulation frequency, the spatial distribution of active cells is generally preserved.Blurring and coarse patterns of activity are observed.Fig. 11 illustrates that at the frequency of 5 Hz, the amplitude of the oscillations is sharply reduced, but the activity region is the largest.

Figure 1 :Figure 2 :
Figure 1: Building an ROI and calibrating the geometric size of pixels

Figure 3 :Figure 4 :
Figure 3: Determination of the radius (R) and angle (Θ): O is the geometric center; C is the center of mass of brightness.Unit of x and y is pixel

Figure 5 :Figure 6 :
Figure 5: The shape of the pulses of mechanical cell contractions when stimulating oscillations of 1 Hz and 4.5 V are applied, estimated for different values of β

Figure 7 :
Figure 7: Change in the rate of contraction of cardiomyocytes in the microchip, depending on the voltage and excitation frequency for the control group (A).R 2 is the square of the correlation coefficient between the approximation function and the experimental data

Figure 8 :
Figure 8: Example of a wave process in a microchip (Condition B, stimulation 1 Hz 4 V).The wave propagates along the X axis, which leads to a delay in cell contraction in the strata.The distance between strata is 380 μm, the width of strata is 54 μm.The blue curve is the oscillations in left strata, the red curve is the oscillations in the right strata

Figure 9 :Figure 10 :
Figure 9: The spread of the average rate of cell contraction in the strata along the axis versus time: (a) along X axis (bright strips that correspond to the maximum rate of the beat are inclined); (b) along Y axis