COMSOL 6.4 - Cell Membrane Electroporation

Electroporation is a technique that employs localized electric fields to generate pores in cell membranes via the rearrangement of the lipid bilayer, improving the cell permeability for ions and pharmaceuticals. High-intensity, nanosecond-short, electric pulses have been shown capable of generating transient nanometric pores onto cell membranes with minimal side effects for the cell. A nanosecond short pulse carries a spectral content in the GHz regime of the frequency domain. Around this frequency, cell membranes are subjected to another concurrent relevant phenomenon called dielectric dispersion, which has been shown to contribute to electroporation efficacy. Both effects are accounted for in this tutorial.

This Cell Membrane Electroporation tutorial model shows the electroporation phenomenon on a standard spherical cell both in the Time Domain and in the Frequency Domain. The model uses an analytical formulation for the electroporation contribution whose parameter values are taken from Ref. 1. In particular, the electroporation effect is accounted for by considering a time-varying conductivity for the membrane:

(1)

where σm0 is the static baseline membrane conductivity, N(t, Vm(t)) is the time-dependent and transmembrane voltage-dependent pore density generation, σp is the pore conductivity, rp is the pore radius, and K is a term proportional to the entrance length and energy barrier of the pore.

On the other hand, the dielectric dispersion is included in this tutorial model by means of a multiple Debye poles dispersion model:

(2)

where εr(f) is the frequency-dependent relative permittivity, ε∞ is the relative permittivity at the high frequency limit, Δεm is the mth relative permittivity contribution, τvm is the relaxation time corresponding to the mth relative permittivity contribution, and N is the number of poles. These dispersive material quantities are taken from Ref. 2 and Ref. 3.

Results such as the electric field, transmembrane voltage, pore density, and membrane conductivity are investigated. They are comparable with the results shown in Ref. 4.

Model Definition

Figure 1 shows the geometry of the cell membrane. The model is solved in 2D axisymmetry, and consists of a spherical cell (that is, a thin insulating spherical membrane layer confining a conductive intracellular domain) placed inside a conductive extracellular domain cylinder. The excitation for electroporation is exerted from the top face of the cylinder whilst the bottom face is grounded.

Figure 2 shows the mesh. It is built using a distribution of elements along the membrane’s boundaries for a lighter, yet accurate, mesh in the proximity of the thin membrane domain. This saves computational costs without penalizing the accuracy of the simulation.

The implementation requires the Electric Currents interface and the Boundary ODEs and DAEs interface. The former interface models the dielectric dispersion of the materials using a multiple Debye poles dispersion subfeature. This effect is shown in Figure 3 by the red dashed line for the membrane, and by the blue solid/green dotted lines for the extra/intracellular domains. The latter interface allows to incorporate on the cell membrane the electroporation-induced conductivity increase (see the Results and Discussion section for details).

Figure 1: 2D axisymmetric geometry of the spherical cell surrounded by the extracellular domain. The cell membrane domain is not visible given its small thickness.

Figure 2: Computational mesh used in this tutorial.

Figure 3: Electrical conductivity (left axis) and relative permittivity (right axis) spectra of all materials used in the model.

Results and Discussion

Two studies are performed in this tutorial: a Small-Signal Analysis, Frequency Domain Perturbation study and a Time Dependent study. Both studies include a first Stationary study step to initialize all variables. Figure 4 and Figure 5 show the frequency response of the system to a 10 V AC perturbation between 1 kHz and 1 THz. Figure 6–Figure 10 show the transient response of the system to a 650 V Gaussian pulse centered at t = 5 ns.

Figure 4 shows the 2D axisymmetric electric potential surface plots and streamlines at three different frequencies. At low frequency (1000 kHz, left inset) the membrane behaves like as a DC capacitor and a relatively small excitation signal couples to the inside of the cell. At large frequency (1 GHz, right inset) the membrane behaves like as an AC capacitor and the excitation signal strongly couples to the intracellular domain. The central inset shows an intermediate condition between small and large frequencies.

Figure 5 displays the transmembrane voltage modulus and phase spectra at a point located at the top of the cell.

Figure 6 depicts the same quantities as Figure 4 but for three time instants. It confirms that a nanosecond short pulse reproduces the qualitative behavior of a high-frequency AC excitation (Figure 4, right inset).

Figure 7 shows the transmembrane voltage, the membrane conductivity, and the pore density at a point located at the top of the cell. The increase of the transmembrane voltage due to the charging of the capacitive membrane triggers an exponential increase in the pore density. This translates into a significant increase in the membrane conductivity that in turn creates a conductive path for the potential built onto the external membrane boundary toward the intracellular domain. As a result, the transmembrane potential drops.

Figure 8 shows the membrane conductivity along the cell’s edge at three time instants. Before 5 ns, the membrane conductivity remains close to a small static value. After 5 ns, the membrane conductivity at the cell’s poles increase significantly and remain to high values for a few more nanoseconds.

Figure 9 is similar to Figure 8 but it displays the transmembrane voltage. After 5 ns, the transmembrane voltage at the cell’s poles drops due to the locally increased membrane conductivity.

Figure 10 shows the 2D-revolved view of the membrane conductivity at the peak vale of the pulse, remarking that the maximum increase of conductivity occurs at the cell’s poles.

Figure 4: 2D electric potential surface plots and streamlines at different frequencies.

Figure 5: Left axis: transmembrane voltage modulus. Right axis: phase spectra.

Figure 6: 2D electric field surface plot and streamlines at different time instants.

Figure 7: Left axis: transmembrane voltage and membrane conductivity (magnified 50 times for better visualization). Right axis: pore density.

Figure 8: Membrane conductivity profile along the cell’s edge at different time instants.

Figure 9: Transmembrane voltage profile along the cell’s edge at different time instants.

Figure 10: 2D-revolved view of the membrane conductivity at the peak value of the pulse.

References

1. G. Pucihar, D. Miklavcic, and T. Kotnik, “A Time-Dependent Numerical Model of Transmembrane Voltage Inducement and Electroporation of Irregularly Shaped Cells,” in IEEE Transactions on Biomedical Engineering, vol. 56, no. 5, pp. 1491–1501, 2009.

2. R. Büchner, G.T. Hefter, and P.M. May, “Dielectric Relaxation of Aqueous NaCl Solutions,” J. Phys. Chem. A, vol. 103, no. 1, pp. 1–9, 1999.

3. B. Klösgen, C. Reichle, S. Kohlsmann, and K.D. Kramer, “Dielectric Spectroscopy as a Sensor of Membrane Headgroup Mobility and Hydration,” J. Biophys., vol. 71, no. 6, pp. 3251–3260, 1996.

4. E. Salimi, Nanosecond Pulse Electroporation of Biological Cells: The Effect of Membrane Dielectric Relaxation, master’s thesis, Dept. Electrical and Computer Eng., Univ. Manitoba, Winnipeg, 2011.

The parameter values of this tutorial model were taken from Ref. 1 with the permission of the authors’ research group (Laboratory of Biocybernetics, Department of Biomedical Engineering, Faculty of Electrical Engineering, University of Ljubljana, Slovenia), as well as from Ref. 2 and Ref. 3. The results of the present model are in agreement with those shown in Ref. 4.

ACDC_Module/Devices,_Resistive/cell_membrane_electroporation

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select > > .

Click .

In the tree, select > > .

Click .

Click

In the tree, select > > .

Click

First, define the geometric and physical parameters of the model.

Global Definitions

Geometric Parameters

In the window, under click .

In the window for , type Geometric Parameters in the text field.

Locate the section. In the table, enter the following settings:


Name	Expression	Value	Description
r_cell	10[um]	1E-5 m	Radius of the cell
t_m	5[nm]	5E-9 m	Thickness of the cell membrane
r_p	0.8[nm]	8E-10 m	Pore radius
h_dom	100[um]	1E-4 m	Height of the extracellular domain
r_dom	3*r_cell	3E-5 m	Radius of the extracellular domain

Physical Parameters

In the toolbar, click

and choose > .

In the window for , type Physical Parameters in the text field.

Locate the section. In the table, enter the following settings:


Name	Expression	Value	Description
Vep	170[mV]	0.17 V	Characteristic voltage of electroporation
sigma_p	0.16[S/m]	0.16 S/m	Electrical conductivity inside the pore
N0	1.5e9[m^-2]	1.5E9 1/m²	Equilibrium pore density
n	0.15	0.15	Relative entrance length of pores
q	2.46	2.46	Pore creation rate
alpha	1e9[m^-2/s]	1E9 1/(m²·s)	Pore creation rate density
w0	2.65	2.65	Pore energy barrier
sigma_m	5e-7[S/m]	5E-7 S/m	Electrical conductivity of the membrane (static value)
Ginf_m	13.9e-12[F/m]/epsilon0_const	1.5699	Relative electrical permittivity at high frequency of the membrane
Gvm1_m	2.3e-11[F/m]/epsilon0_const	2.5976	First relative electrical permittivity contribution of the membrane
Gvm2_m	7.4e-12[F/m]/epsilon0_const	0.83576	Second relative electrical permittivity contribution of the membrane
tauvm1_m	3e-9[s]	3E-9 s	First relaxation time of the membrane
tauvm2_m	4.6e-10[s]	4.6E-10 s	Second relaxation time of the membrane
sigma_e	0.14[S/m]	0.14 S/m	Electrical conductivity of the extracellular domain
sigma_i	0.3[S/m]	0.3 S/m	Electrical conductivity of the intracellular domain
Ginf_f	5.9e-10[F/m]/epsilon0_const	66.635	Relative electrical permittivity at high frequency of cellular domains
Gvm1_f	1.18e-10[F/m]/epsilon0_const	13.327	Relative electrical permittivity contribution of cellular domains
tauvm1_f	6.2e-12[s]	6.2E-12 s	Relaxation time of the cellular domains

Study Parameters

In the toolbar, click

and choose > .

In the window for , type Study Parameters in the text field.

Locate the section. In the table, enter the following settings:


Name	Expression	Value	Description
V_fd	10[V]	10 V	Amplitude of the frequency domain perturbation
V_td	650[V]	650 V	Amplitude of the time dependent pulse
t_end	10[ns]	1E-8 s	Duration of the time dependent study

Build the geometry following these steps.

Geometry 1

Cell

In the window, expand the > node.

Right-click and choose .

In the window for , locate the section.

In the text field, type r_cell.

In the text field, type 180.

Locate the section. In the text field, type -90.

In the text field, type Cell.

Click to expand the section. In the table, enter the following settings:


Layer name	Thickness (m)
Layer 1	t_m

Cellular Domains

In the toolbar, click

In the window for , type Cellular Domains in the text field.

Locate the section. In the text field, type r_dom.

In the text field, type h_dom.

Locate the section. In the text field, type -h_dom/2.

Form Union (fin)

In the window, click .

In the window for , click

Click the

button in the toolbar.

Define useful selections under the Geometry node.

Membrane

In the toolbar, click

and choose .

In the window for , type Membrane in the text field.

On the object , select Domains 2 and 4 only.

Intracellular Domain

In the toolbar, click

and choose .

On the object , select Domain 3 only.

In the window for , type Intracellular Domain in the text field.

Click

Extracellular Domain

In the toolbar, click

and choose .

On the object , select Domain 1 only.

In the window for , type Extracellular Domain in the text field.

Intra/Extracellular Domains

In the toolbar, click

and choose .

On the object , select Domains 1 and 3 only.

In the window for , type Intra/Extracellular Domains in the text field.

Membrane Boundary Intra

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

In the text field, type Membrane Boundary Intra.

On the object , select Boundaries 12 and 13 only.

Membrane Boundary Extra

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

In the text field, type Membrane Boundary Extra.

On the object , select Boundaries 11 and 14 only.

Terminal Boundary

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

On the object , select Boundary 8 only.

In the text field, type Terminal Boundary.

Ground Boundary

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

In the text field, type Ground Boundary.

On the object , select Boundary 2 only.

Set up two general extrusions to map the potential between membrane boundaries; use them to define the membrane potential difference.

Definitions

Mapping from Intracellular to Extracellular

In the window, expand the > node.

Right-click and choose > .

In the window for , type Mapping from Intracellular to Extracellular in the text field.

In the text field, type in2out.

Locate the section. From the list, choose .

From the list, choose .

Locate the section. In the text field, type r/sqrt(r^2+z^2)*(r_cell-t_m).

In the text field, type z/sqrt(r^2+z^2)*(r_cell-t_m).

Locate the section. Select the checkbox.

Mapping from Extracellular to Intracellular

In the toolbar, click

and choose .

In the window for , type Mapping from Extracellular to Intracellular in the text field.

In the text field, type out2in.

Locate the section. From the list, choose .

From the list, choose .

Locate the section. In the text field, type r/sqrt(r^2+z^2)*r_cell.

In the text field, type z/sqrt(r^2+z^2)*r_cell.

Locate the section. Select the checkbox.

Transmembrane Potential

Right-click and choose .

In the window for , type Transmembrane Potential in the text field.

Locate the section. From the list, choose .

From the list, choose .

Locate the section. In the table, enter the following settings:


Name	Expression	Unit	Description
Vm	out2in(V)-V	V

Membrane Variables

Right-click and choose .

In the window for , locate the section.

From the list, choose .

In the text field, type Membrane Variables.

Locate the section. In the table, enter the following settings:


Name	Expression	Unit	Description
vm	in2out(Vm)*F_const/R_const/minput.T		Dimensionless transmembrane voltage
K	((exp(vm)-1)/(exp(vm)(w0exp(w0-nvm)-nvm)/(eps+w0-nvm) - (w0exp(w0+nvm)+nvm)/(w0+n*vm)))		Electroporation parameter related to wo and n
sigma_m_t	in2out(N)sigma_ppir_p^2K		Time-dependent membrane conductivity

Gaussian Pulse

In the toolbar, click

and choose .

In the window for , type Gaussian Pulse in the text field.

Locate the section. In the text field, type t_end/2.

In the text field, type t_end/25.

From the list, choose .

In the text field, type V_td.

Use the following settings for the Electric Current interface and the Boundary ODE interface.

Electric Currents (ec)

Current Conservation - Membrane

In the window, under > click .

In the window for , type Current Conservation - Membrane in the text field.

Locate the section. From the list, choose .

Dispersion 1

In the window, click .

In the window for , locate the section.

From the list, choose .

In the N text field, type 2.

Current Conservation - Intra/Extracellular

In the toolbar, click

and choose .

In the window for , type Current Conservation - Intra/Extracellular in the text field.

Locate the section. From the list, choose .

Dispersion 1

In the window, click .

In the window for , locate the section.

From the list, choose .

In the N text field, type 1.

Ground 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Terminal for Small-Signal Analysis, Frequency Domain

In the toolbar, click

and choose .

In the window for , type Terminal for Small-Signal Analysis, Frequency Domain in the text field.

Locate the section. From the list, choose .

Harmonic Perturbation 1

In the toolbar, click

and choose .

In the window for , locate the section.

In the V0 text field, type V_fd.

Terminal for Stationary Before Time Dependent

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. From the list, choose .

In the V0 text field, type 1e-9.

In the text field, type Terminal for Stationary Before Time Dependent.

Terminal for Time Dependent

In the toolbar, click

and choose .

In the window for , type Terminal for Time Dependent in the text field.

Locate the section. From the list, choose .

In the V0 text field, type gp1(t).

Boundary ODEs and DAEs (bode)

In the window, under click .

In the window for , locate the section.

From the list, choose .

Locate the section. Click

In the table, enter the following settings:


Dependent variable quantity	Unit
Custom unit	1/m^2

In the table, enter the following settings:


Source term quantity	Unit
Custom unit	1/m^2/s

Click to expand the section. In the text field, type N.

In the table, enter the following settings:


N

Distributed ODE 1

In the window, under > click .

In the window for , locate the section.

In the f text field, type alpha*exp((Vm/Vep)^2)*(1-N/N0*exp(-q*(Vm/Vep)^2)).

Initial Values 1

In the window, click .

In the window for , locate the section.

In the N text field, type N0.

Add one blank material per domain and link the corresponding dielectric dispersion properties.

Materials

Cell Membrane

In the window, under right-click and choose .

In the window for , locate the section.

From the list, choose .

In the text field, type Cell Membrane.

Locate the section. In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Electric conductivity	sigma_iso ; sigmaii = sigma_iso, sigmaij = 0	sigma_m+sigma_m_t	S/m	Basic
Relative permittivity	epsilonr_iso ; epsilonrii = epsilonr_iso, epsilonrij = 0	Ginf_m+Gvm1_m+Gvm2_m	1	Basic
Relaxation times	{tauvm1, tauvm2}	{tauvm1_m, tauvm2_m}	s	Multipole Debye
Relative permittivity contributions	{Gvm1, Gvm2}	{Gvm1_m, Gvm2_m}	1	Multipole Debye

Intracellular Domain

Right-click and choose .

In the window for , locate the section.

From the list, choose .

In the text field, type Intracellular Domain.

Locate the section. In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Electric conductivity	sigma_iso ; sigmaii = sigma_iso, sigmaij = 0	sigma_i	S/m	Basic
Relative permittivity	epsilonr_iso ; epsilonrii = epsilonr_iso, epsilonrij = 0	Ginf_f+Gvm1_f	1	Basic
Relaxation times	tauvm	tauvm1_f	s	Multipole Debye
Relative permittivity contributions	Gvm	Gvm1_f	1	Multipole Debye

Extracellular Domain

Right-click and choose .

In the window for , locate the section.

From the list, choose .

In the text field, type Extracellular Domain.

Locate the section. In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Electric conductivity	sigma_iso ; sigmaii = sigma_iso, sigmaij = 0	sigma_e	S/m	Basic
Relative permittivity	epsilonr_iso ; epsilonrii = epsilonr_iso, epsilonrij = 0	Ginf_f+Gvm1_f	1	Basic
Relaxation times	tauvm	tauvm1_f	s	Multipole Debye
Relative permittivity contributions	Gvm	Gvm1_f	1	Multipole Debye