Simulation of single cell electroporation

. This paper attempts to introduce the dynamics of electroporation into the single cell model. The main characteristics of the model were described. Results show the generation and rise of pores in a round cell with a radius of 50 μm exposed to 40 kV/m electric field for 1 ms were analysed, so that the transmembrane potential, the number of pores and the membrane conductance could be calculated. Finally, how the model can help explain the experiment is discussed.


Introduction
Electroporation is a kind of microbiological technology. With cells subject to electric field pulse, it is easier for osmosis, making chemicals or DNA and other foreign substances enter cells. In an electric field, there comes potential distribution varying from space to space [1]. In order to obtain the plasma membrane permeability, the electric field intensity above the lower threshold is needed to achieve the effective transmembrane potential, so as to open the pores on the plasma membrane. Electroporation, especially reversible electroporation, is a complex and costly project. If the electric field intensity is higher than the upper limit of the threshold, it will lead to irreversible cell infiltration or even death [2], which is often used in some application scenarios like tumour resection. Sung et al. [3] suggested that the clinical results of the irreversible electroporation on a liver tissue could be predicted through mathematical calculation.
It is necessary to supplement experimental knowledge with theoretical models. Mathematical modelling is a method to study the electrical phenomenon in the process of electroporation, which can speed up the research process and save the cost. The calculation of electric field distribution is a relatively simple and effective tool to analyse and interpret the experimental results.
There are two kinds of electroporation models. The first model [4][5][6][7] focuses on reproducing the time course of pore generation and evolution, while it simplified the condition of spatial constraints. The second model [8][9][10][11][12][13][14][15] takes the opposite approach, valuing potential field boundary issue more than temporal consideration. However, in these studies, the timing of electroporation is either very simple or not mentioned at all. This paper attempts to introduce the dynamics of electroporation into the single cell model which equally focuses on time and space, and can study electroporation in millisecond time scale. The main characteristics of the model were described. Results show the generation and rise of pores in a round cell with a radius of 50 μm exposed to 40 kV/m electric field for 1 ms were analysed, so that the transmembrane potential, the number of pores and the membrane conductance could be calculated. Finally, how the model can help explain the experiment is discussed.

Numerical simulation model
In this study, a round cell with radius a was exposed to an electric field E . We define two indices, e  and i  , which stand for intracellular and extracellular potentials .
where  being the distance from the cell center,  being the polar angle measured relative to the electric field E direction. For the current density, in which n being the unit vector perpendicular to membrane surface. m V is the transmembrane potential, which is obtained by subtracting the transmembrane potential from the intramembrane potential. s stands for conductivity. 2 2 10 F/m m C   being the surface capacitance of the membrane [15]. Since the change of m C measured in the experiment is 2% [16], it is assumed to be constant. 2 2 S/m g   being the surface conductance of the membrane [15]. rest 0.08 V V   being the rest potential [17].
In the simulation, the cell surface is discretized from 0 to  in steps of   in  . A  indicates a small fraction of cell area covered by   . Supposing there are K  pores in A  , the current of all pores is regarded as parallel connection, and the total current in A  is obtained by accumulation. Therefore, we have where p R is the resistance across the pore and i R is the input resistance, both of which are considered to be in series.
In the initial stage of 0 t  and 0 E  , the cell is at its resting potential rest V , and the extracellular potential is 0, which is the initial condition. For 0 t  , E is set non-zero and the potential evolves over time according to equation (3). N . Therefore, the pore density begins to decrease if the pulse is shut down.

Numerical Simulation
We built a model as shown in Fig. 1A. And fig. 1B shows the simulated electric field environment. The spherical coordinates are used to discretize the intracellular and extracellular space. The discrete step size is / 128 (1) is transformed into a set of linear equations by finite difference method. As a Dirichlet boundary condition, the external electric field is applied to the potential e  along the outer boundary of the cell shell layer ( 3a   ). According to the continuity of the membrane, the current value is set to p I  in equation (3). For discretized elements, the superscript is related to quantities that vary with the polar angle. The number of pores is calculated for each discrete angle on the membrane. The resulting p I  will be used to determine intracellular and extracellular potentials in the next time step.
Five locations on the cell membrane are particularly important: depolarization pole D , hyperpolarization pole H , equator E , and the boundary between electroporation and non-electroporation regions ( b D and b H ). Fig. 1C and Fig. 1D show the electric field and potential distribution of the simulated cells, respectively.  Fig. 2 shows a 50 μm cell exposed to an electric field of 40 kV/m for 1 millisecond. The charging phase continued until the pore forming stage started. When the membrane does not yet contain any pores, the membrane charge is similar to that of a parallel resistor and capacitor. In Fig. 2A, within the initial 0.51 μs, / rest The voltage amplitudes at these locations also decrease considerably. It is found that b D m V decreases at 1.5 μs, and the pore will not form in this area immediately, which indicates that the transmembrane potential will be affected by the state of the adjacent area, resulting in the decrease of transmembrane potential, even in the area of cell membrane without electroporation. When the pore growth rate is less than 10 -6 per time step, the pore nucleation stage is completed. Taking   At 1 ms, the transmembrane potential was stable.

Results
Electroporation reduces the amplitude of m V  , not only in the area of local pores, but also in the whole cell range. The interaction of pore number and pore size makes an interesting change in the conductance of cell membrane (see Fig. 3C). In the early stage of electroporation, the surface conductivity G  is the largest. The increase of G  is more significant in the area where lots of pores come into being, mainly concentrated near the two poles instead of the boundary area. With the evolution of pore radius, the amplitude of G  distribution decreases and becomes wider.

Discussion
Because the cell parameter settings in different experiments often come from different literature results, the performance of the model will show some differences, which makes it difficult to reproduce the results of the paper. We chose settings from Hibino et al. [15] to match the parameters. Other parameters came from other sources, such as experiments about artificial lipid bilayers. In fact, the mathematical model cannot replace the experimental work, it is only an additional source of information to explain the phenomenon and the experimental plan. In a word, it is helpful to understand the process of electroporation.
The most detailed information on the process of electroporation in cells comes from the transmembrane potential measurements obtained by Hibino et al. [15] using voltage sensitive fluorescent dyes. By observing the changes of fluorescence around cells during electroporation, we can understand its spatial distribution and temporal evolution. The model reproduces the transmembrane potential well. As we can see in the experiment, the model fluctuates in microsecond time scale, and the transmembrane potential shows the saturation characteristic of electroporation. In the model and experiment, due to the negative electrostatic potential, the formation of the pores first occurs at the hyperpolarized pole. There is no such slow depolarization drift of cell potential in the model. As we can see, after the first three milliseconds, the poles change very little.
Some studies used equation (6) to charge passivated cells to assess which parts of the cells were electroporated [9,10,22,23]. The potential hypothesis is that the membrane will be electroporated regardless of the transmembrane potential found in equation (6). The generation of pores reduces the number of pores anywhere on the membrane, not just in the electroporation area. The cells never reach the transmembrane potential predicted by equation (6), so the degree of electroporation is smaller.
Compared with the conductivity of cytoplasm and extracellular medium, the conductivity of plasma membrane is very low. In this case, the conductivity of plasma membrane increases, which makes the distribution of electric field more complex. Because of the stomata, the membrane conductance increases sharply, which makes the cell unable to maintain its negative resting potential, so the equatorial transmembrane potential is very small. In order to keep the balance of current input and output, the conductance of the two hemispheres should be approximately equal. However, Hibino et al. [15] inferred that the membrane conductance is asymmetric: D G is about twice that of H G . This asymmetry may be due to the approach of transmembrane voltage to positive value, which leads to cell depolarization. As mentioned above, the model does not reproduce this drift.
The electroporation stage observed in the model is very similar to that determined by Teissie et al. [23]. Although some predictions of the model still need to be confirmed by experiments, the similarity between the electroporation models observed in the experiments convinces us that the single cell model proposed in this paper may facilitate practical application.
Vitro electroporation usually involves only one type of tissue. It is easy to produce non-uniform electric field because of the flat electrode placed on the surface of tissue or needle electrode inserted into tissue, and the shape of tissue is also different. In addition, biological tissue composed of cells cannot be regarded as a uniform conductor. Compared with the conductivity of cytoplasm and extracellular medium, the conductivity of plasma membrane is very low. In the case of electroporation, the conductivity of plasma membrane increases, which greatly changes the electrical characteristics in tissues and makes the distribution of electric field more complex.
In the case of electroporation in vivo, different tissues and organs have different electrical properties. In addition, the diversity of their structures leads to the high inhomogeneity of the electric field generated by pulse transmission. The diversity of the shapes, sizes and positions of electrodes used for electroporation in vivo and the resulting differences in biological reactions further indicate that it is necessary to study the differences in electric field distribution caused by specific electroporation.

Conclusion
The purpose of this paper was to figure out a model for single cell electroporation. Although some parameters of the model still need to be further confirmed, we are convinced that these results supported the potential practical applications of this single cell model proposed in this paper. This makes it possible to facilitate electroporation protocols for specific tasks, such as creating pores for corresponding functions like DNA uptake or drug delivery.