OpenEP: an open-source simulator for electroporation-based tumor treatments

Electroporation (EP), the increase of cell membrane permeability due to the application of electric pulses, is a universal phenomenon with a broad range of applications. In medicine, some of the foremost EP-based tumor treatments are electrochemotherapy (ECT), irreversible electroporation, and gene electrotransfer (GET). The electroporation phenomenon is explained as the formation of cell membrane pores when a transmembrane cell voltage reaches a threshold value. Predicting the outcome of an EP-based tumor treatment consists of finding the electric field distribution with an electric threshold value covering the tumor (electroporated tissue). Threshold and electroporated tissue are also a function of the number of pulses, constituting a complex phenomenon requiring mathematical modeling. We present OpenEP, an open-source specific purpose simulator for EP-based tumor treatments, modeling among other variables, threshold, and electroporated tissue variations in time. Distributed under a free/libre user license, OpenEP allows the customization of tissue type; electrode geometry and material; pulse type, intensity, length, and frequency. OpenEP facilitates the prediction of an optimal EP-based protocol, such as ECT or GET, defined as the critical pulse dosage yielding maximum electroporated tissue with minimal damage. OpenEP displays a highly efficient shared memory implementation by taking advantage of parallel resources; this permits a rapid prediction of optimal EP-based treatment efficiency by pulse number tuning.


Numerical model
In what follows some details of the numerical model discussed in the Methods section are presented.
The non-linear Laplace's equation (Eq. 1) has three terms: The first term is discretized with the following central difference scheme: Then, each sub-term is again discretized with a central scheme: σ i,j,k + σ i−1,j,k 2 Φ i,j,k − Φ i−1,j,k ∆x 1 ∆x The y and z terms of the Eq. A are treated analogously. The resulting equation is solved for Φ i,j,k through a standard relaxation procedure. The electric field is computed through the gradient of the electrostatic potential with a central difference scheme: where the derivative of Φ is approximated as: The bioheat equation is written: The first three terms from the right side of the Eq. D are similar to those in Eq. A, therefore having identical treatment, for instance, the first term becomes: the term σ |∇Φ| 2 in the Eq. D is approximated in a similar way to that in Eq. C. The final discretization reads: The previous finite difference scheme was solved for (T n+1 i,j,k ) with an standard explicit method. To compute the electric current (Eq. 6), the domain is partitioned by half, leaving one electrode in each half-domain. The integral is computed across the surface of one of the halves. For example, to calculate the outgoing flow through one of the half-domain surfaces (when x = ii/2), the discretization is computed as: Analogously, to calculate the outgoing flow when x = 0: These calculations are similar for y and z axis.

Boundary Conditions
For every electrode node at the upper surface (k = 0), Eq. 8 was discretized using a first order scheme: The other boundaries are six rectangular surfaces exposed to the air (due to the box-shape domain). Temperature and electric potential are subject to the Neumann boundary conditions according to Eq. 9 and Eq. 10. Particularly, for the left surface, i.e. an yz plane, the boundary conditions are: Recalling that κ = 0 and using a first order scheme in the approximation of the derivatives of Eq. G, the following results: T 0,j,k = T 1,j,k . The remaining boundaries are computed analogously.
Compilation and running script options NOTE: All files *.cpp and *.h used to generate the simulation results can be accessed at https://github.com/LSC-UBA/OpenEP/tree/master/src run.sh options Description -d <simulation-directory> Define the simulation directory. Default: simulation-i, where i is the next number regarding to the previous simulation directory.

-n <number-of-threads>
Define the number of threads that will be used by the simulation. Default: number of logical cores of the computer node.

-c <compilation-option>
Define the compilation option (see Table S.2). Default: FAST_OMP.   ("vtk" or "csv") save_format "vtk" "vtk" How to optimize an EP-based treatment in terms of pulse number in five steps.    Having the electroporated tissue trajectory, it is possible to determine the pulse dosage maximizing electroporated tissue with minimum damage (if the damage is not included, as in the green curve, the optimum pulse dosage (pulse number) is obtained as the first pulse in which the maximum difference between EA (n+1) -EA(n) is smaller than a prescribed small value). It is optimal because further pulsing will not increase the electroporated tissue area but may introduce damage.