A fractional stochastic theory for interfacial polarization of cell aggregates
AA fractional stochastic theory for interfacial polarization of cell aggregates
Pouria A. Mistani , ∗ Samira Pakravan , and Frederic Gibou , , Department of Mechanical Engineering, Department of Computer Science, Department of Mathematics,University of California Santa Barbara,Santa Barbara, CA 93106, USA (Dated: August 28, 2020)We present a theoretical framework to model the electric response of cell aggregates. We establisha coarse representation for each cell as a combination of membrane and cytoplasm dipole moments.Then we compute the effective conductivity of the resulting system, and thereafter derive a Fokker-Planck partial differential equation that captures the time-dependent evolution of the distributionof induced cellular polarizations in an ensemble of cells. Our model predicts that the polarizationdensity parallel to an applied pulse follows a skewed t-distribution, while the transverse polariza-tion density follows a symmetric t-distribution, which are in accordance with our direct numericalsimulations. Furthermore, we report a reduced order model described by a coupled pair of ordinarydifferential equations that reproduces the average and the variance of induced dipole moments inthe aggregate. We extend our proposed formulation by considering fractional order time derivativesthat we find necessary to explain anomalous relaxation phenomena observed in experiments as wellas our direct numerical simulations. Owing to its time-domain formulation, our framework can beeasily used to consider nonlinear membrane effects or intercellular couplings that arise in severalscientific, medical and technological applications.
I. INTRODUCTION
Effects of external electric fields on heterogeneous sys-tems have been of great scientific and technological im-portance throughout the past century. These systemsinclude composite materials, colloidal suspensions, andbiological cells. In the case of biological cells, electricfields have found several applications for cell fusion [1, 2],electrorotation [3, 4], dielectrophoresis [5, 6] and cancercell separation [7, 8], electroporation [9], levitation [10]and cell deformation [11]. For an early review of its appli-cations in biotechnology and medicine we refer to Markxand Davey (1999) [12]. More recently, transmembranepotential (TMP) patterns that emerge in multicellularliving organisms have gained extra attention due to dis-covery of their regulatory role in development and regen-eration; we refer to the review of Levin et al. (2017) for acomprehensive overview on developmental bioelectricity[13]. Even though modern research on TMP manipula-tions are focused on molecular based treatments, it hasbeen long known that TMP patterns altered by externalelectric fields could influence development [14], embryo-genesis [15, 16] or wound healing [17]. Therefore, devel-oping a predictive and generalizable mathematical modelto understand the effects of different cellular mechanismson the tissue level bioelectric patterns poses promisingopportunities in bioengineering.In all of these applications, it is essential to have a pre-cise knowledge of the TMP induced over cell membranes,especially in cell aggregates that are composed of tens ∗ Corresponding author: [email protected] of thousands of cells with a heterogeneous mix of mor-phologies and electrical properties. Even though muchis known about single cell electric interactions, a theorythat predicts a detailed distribution of transmembranepotentials within cell aggregates has been missing. In thiswork, we present a novel theoretical framework based ona Fokker-Planck description, for tracing the time evolu-tion of the probability density of multicellular polariza-tions in response to arbitrary electric stimulations. Wefurther introduce a moment-based analytic reduced-ordermodel of the proposed Fokker-Planck equation that pro-vides statistics of transmembrane potentials with mini-mal computational expense. Importantly, besides multi-cellular systems, our theory is applicable to a broad rangeof systems such as shelled colloidal particles, emulsionsand composite materials [18]. Moreover, it can be ex-tended to include the effects of membrane nonlinearities[19] or intercellular couplings [20], counterion polariza-tions, and eventually real-time pulse optimization, whichis of great benefit to emerging biomedical treatments us-ing electric fields such as electrochemotherapy.
A. Physical bioelectric processes
Electrical properties of biological tissues have beenextensively investigated for the last two centuries sincethe discovery of Ohm’s law. A comprehensive histori-cal overview of different aspects of biological dielectricresponse (including basic concepts, tabulated data, un-derlying molecular processes and effective cellular in-teractions) is covered in the surveys of Schwan (1957)[21], Stuchly and Stuchly (1980) [22], Pethig (1984) [23],Pethig and Kell (1987) [24], Foster and Schwan (1989) a r X i v : . [ c s . C E ] A ug [25], McAdams and Jossinet (1995) [26], Gabriel (1996)[27], Kuang and Nelson (1998) [28], with more recent re-views on its different applications such as electroporationprovided by Kotnik et al. (2019) [29].Early research on the electric response of bulk biolog-ical materials revealed that tissues exhibit resistive andcapacitive behaviors. Experiments showed an early peakin current in response to a step voltage, which couldbe attributed to an increasing tissue resistance, or froman induced counter-potential polarization. Furthermore,with the advent of high frequency apparatus it was pos-sible to examine the high frequency response of tissues,which led to the recognition that tissues exhibit low re-sistance at higher frequencies and the dielectric responseof tissues was determined by different physical processesat each frequency regime. Bulk electric properties aremainly determined by cell membranes and cellular struc-tures. At the cellular level, three main physical processes,namely interfacial polarization, ionic diffusion and dipo-lar orientation of polar molecules, are identified to playkey roles in dielectric dispersions at different frequencyregimes [21]. While the origin of α - and γ -dispersions arerelaxation processes in the bulk phases of the material, β -dispersions originate from internal boundary conditionsimposed by interfaces separating different phases; this isthe focus of the current work. Below, we briefly reviewthese mechanisms. • α -dispersion: The main factor that contributes tothe α -dispersion (at audio or sub-KHz frequencies) is theionic diffusion in the electric double layer in the immedi-ate vicinity of charged surfaces as well as in the bulk. α -dispersion is characterized with high dielectric constantsat low frequencies. Schwan first observed this mode ofpolarization at low frequencies in biological tissues [21]and later, Schwan and co-workers showed that this effectis also observed in non-biological colloids [30]. Schwarz(1962) [31] was the first to develop a theory that tookinto account the counterion polarization around colloidalparticles suspended within electrolytes. Schwarz showedthese displacements could be modeled by an additional“apparent” dielectric constant that reach high values atlow frequencies. Schwarz’s method did not consider dif-fusion in the double layer itself, and was later extendedby Einolf and Carstensen (1971) [32, 33] to include diffu-sion on both sides of the membranes. Dukhin and Shilov(1974) [34] proposed a more accurate treatment by con-sidering ionic diffusion in the bulk, rather than just thethin layers around charged particles (see also the reviewby Mandel and Odijk [35]). A simplified formulation ofDukhin’s model that admits analytical solution was givenby Grosse and Foster [36], which helped to show that thecorresponding Cole-Cole plot is broader than the Debye’smodel, indicating that counter-ion polarization is partlyresponsible for the observed anomalous relaxation of bio-logical matter. In short, counter-ion polarization theoriesexplain α -dispersion and predict high permittivities atlow frequencies that exhibit broad relaxation behaviors. • β -dispersion: Interfacial polarization dominates the dielectric properties of tissues at β -dispersion (radiofrequencies from tens of KHz to tens of MHz range,timescales determined by membrane resistance and ca-pacitance) as well as the dielectric properties of colloidsand emulsions. Biological mixtures of interest to ushave a triphasic dielectric structure with conductive partscomposed of a cytoplasm covered by a membrane im-mersed in a continuous medium. Historically, dielectrictheories of interfacial polarization began by considering diphasic suspensions in the seminal treatise of Maxwell(1873) [37] and later Wagner (1914) [38]. In 1925, Fricke[39] developed a dielectric theory for spherical particlessurrounded by nonconductive membranes and extendedit to membrane-covered ellipsoidal particles in 1953 [40].Maxwell-Wagner theory has the following limitations: (i)it is only valid at very low concentrations and assumedthat the local electric field was equal to the externalelectric field, (ii) the interior of particles were assumedto be at constant potential, and (iii) the external fieldwas modeled as if the particle was a perfect insulator.Hanai (1960) [41] developed an interfacial polarizationtheory that was valid at high concentrations, by assum-ing Wagner’s relation holds during successive additionsof infinitesimally small quantities of the disperse phase.Hanai and co-workers (1979) [42] later generalized theirapproach to the case of suspensions of shelled spheres andZhang et al. (1983-1984) [43, 44] showed that the the-ory could explain experiments with polystyrene micro-capsules. This strategy was applied to three-phase struc-tures in 1993 [45]. We shall emphasize that even thoughHanai’s approach is advantageous over Wagner’s theoryas it holds its validity to high concentrations, it is stillbased on the non-conductive assumption for membranesand, more importantly, it is not clear how to considernonlinear variations in the membrane conductance dur-ing the application of electric pulses similar to the caseof electroporation. Our theory builds on this line of workand aims to address these shortcomings by constructinga time-domain model for interfacial polarization in cellaggregates. Unlike the aforementioned works that arelimited to modeling average properties at the aggregatelevel, our approach captures detailed information aboutthe distribution of induced polarizations as well as theirtime-dependent evolution. • γ -dispersion: The third mechanism that is respon-sible for the γ -dispersion (microwave frequencies fromMHz to GHz range) is the dipolar orientations of per-manent polar molecules, e.g. water molecules and othermacromolecules. Under an applied electric torque andopposed by thermal agitations in the medium, polarmolecules undergo rapid reorientations towards thermalequilibrium and exhibit dielectric relaxation. This phe-nomenon is described by Debye’s theory (1929) [46],which is inherently a Fokker-Planck equation that de-scribes the evolution of the probability density of dipo-lar orientations under an applied pulse. Our theory pre-sented here is inspired by this strategy. FIG. 1. Cell membranes are modeled with an array of par-allel resistors and capacitors in a thin layer surrounding thecytoplasm. Here we consider a sharp interface.
B. Equations of interfacial polarization
In its general form, Maxwell’s equations in matter read ∇ · D = ρ f (1) ∇ · B = 0 (2) ∇ × E = − ∂ t B (3) ∇ × H = J f + ∂ t D (4)where E and B are the electric and magnetic fields, re-spectively, D = (cid:15) E , H = µ B and the total current isdefined by J = J f + ∂ t D . However, electric interactionswithin a multicellular system can be modeled under the quasi-electrostatic assumption, i.e. when the wavelengthof the stimulating electric pulse is larger than the cellsize. Under the quasi-electrostatic assumption, the in-duced magnetic fields are negligible and therefore theelectric field is curl free and we may define the electricpotential u by the relation E = −∇ u. Also, comput-ing the divergence of (4) and using (1), we have that ∇ · ( σ ∇ u) = ∂ρ f ∂ t , where in the absence of a net freecharge density ρ f , when only interfacial polarization ispresent, the electric field is given by the solution of theLaplace equation ∇ · ( σ ∇ u) = 0, and we can neglect thepermittivity of the cytoplasm and of the extra-cellularmedium.As for cells, we consider a shelled particle model withtwo concentric surfaces Γ ± forming a membrane withthickness h (see figure 1). The boundary conditions im-pose that the electric potential must be continuous acrossthe interfaces u c ( x − ) = u m ( x − ) , u m ( x + ) = u e ( x + ) , with x ± ∈ Γ ± and that the normal component of the to-tal current J k = ( σ k + j ω(cid:15) (cid:15) k ) E k = Λ ∗ k E k (with k = e , mand c referring to the extra-cellular medium, the mem-brane and the cytoplasm, respectively) must be continu-ous across the boundaries:Λ ∗ c ∂ n u c ( x ) = Λ ∗ m ∂ n u m ( x ) , x ∈ Γ − , Λ ∗ m ∂ n u m ( x ) = Λ ∗ e ∂ n u e ( x ) , x ∈ Γ + . Note that E m · n = − [u] / h. This set of equations havebeen considered by Miles and Robertson (1932) [47] fora single sphere. We further assume a thin membrane bysetting h / R →
0. The exact response of the TMP toa step pulse for an isolated sphere with constant mem-brane conductance has been studied by many authors[48–52]. Even though these equations are not solved forarbitrary cell geometries, Kotnik & Miklav´ci´c [53] provideanalytical solutions for the TMP over oblate and prolatecells in the special case where the cell’s axes of symmetryis parallel to the applied field. Qualitatively, they haveshown that the maximum TMP increases with the equa-torial radius of spheroids, i.e. the radius perpendicularto the applied field. The authors considered an insulat-ing membrane to simplify the analysis; however it wasshown in experiments that this is not a valid assumption[54] and one has to consider changes in membrane pora-tion and permeabilization under an applied pulse. Thisre-structuring of cells membrane under an electric fieldcan be considered, for example, by adopting a nonlinearphenomenological model for the membrane conductance[55]: S m ( t, [ u ]) = S L + S X ( t, [ u ]) + S X ( t, [ u ]) , (5)where S , S and S are the conductance values of themembrane in the resting, porated and permeabilizedstates, respectively. In this model, the level of porationand permeabilization of the membrane are captured inthe functions X and X , which are calculated as a func-tion of the TMP by solving a set of nonlinear ordinarydifferential equations. Overall, we therefore adopt thefollowing boundary value problem to model electric in-teractions in charge-free mixtures, ∇ · ( σ ( x ) ∇ u ) = 0 , x ∈ (Ω c ∪ Ω e ) (6a)with boundary conditions,[ σ ( x ) ∂ n u ] = 0 x ∈ Γ , (6b) C m ∂ t [ u ] + S ( t, [ u ]) [ u ] = σ ( x ) ∂ n u x ∈ Γ , (6c) u ( t, x ) = g ( t, x ) x ∈ ∂ Ω , (6d)and homogeneous initial condition, u (0 , x ) = 0 x ∈ Ω (6e)where we used [ o ] = o e − o c to describe the jump op-erator across the interface Γ in the normal direction, σ c , σ e and σ m are the conductivities of the cytoplasm,the extra-cellular medium and the cell membrane, re-spectively. Note that σ m ≡ S m h and (cid:15) (cid:15) m ≡ C m h aremembrane conductivity and permittivity, respectively. FIG. 2. A snapshot of the 3D spherical tumour composed of ∼ ,
000 random cells considered in this work. Cells are coloredaccording to their TMP values, with redder colors indicating higher positive TMP while bluer colors indicate lower negativeTMP. The Octree data structure is shown on an equatorial slice to emphasize the enhanced resolutions close to cell membranes.
C. Direct numerical simulations
Numerical simulations can be used to directly inves-tigate tissue-level properties of cell aggregates emergingfrom the set of equations 6a–6e along with (5), when con-sidered in a large heterogeneous environment. However,several computational challenges need to be addressedin order to obtain guaranteed accuracy and convergenceof the numerical results, while simultaneously consider-ing large enough number of cells. Particularly, impos-ing jump conditions on numerous irregular sharp inter-faces requires efficient numerical discretization methodsalong with scalable parallel computing algorithms formesh generation and storage as well as advanced linearsystem solvers and preconditioners.In this vein, Guittet et al. (2016) [56] introducedthe Voronoi Interface Method (VIM) to solve Ellipticproblems with discontinuities across the interface of ir-regular domains. Basically, VIM utilizes an interface-fitted Voronoi mesh before applying the dimension-by-dimension Ghost Fluid Method [57]. Importantly, VIMproduces a linear system that is symmetric positive def-inite with only its right-hand-side affected by the jumpconditions. The solution and the solution’s gradients aresecond-order accurate and first-order accurate, respec-tively, in the L ∞ -norm. Later, Guittet et al. (2017) applied VIM to the case of cell electroporation [58] ina serial computing environment. Recently, Mistani etal. (2019) [59] extended these results to parallel comput-ing environments and considered a large tumour spheroidcomposed of ∼ ,
000 ellipsoidal cells. Figure 2 illus-trates a snapshot of this simulation with the transmem-brane potentials depicted over cell membranes. This sim-ulation leveraged the suite of data structures and rou-tines provided by the Portable, Extensible Toolkit forScientific Computation (
PETSc ) [60–62] using the Bi-CGSTAB solver [63] over the linear system precondi-tioned by hypre [64] library. Creation and managementof adaptive octree grids was handled by p4est softwarelibrary [65] along with voro++ [66] library for building anadaptive Voronoi tessellation from the underlying Octreegrid.To our knowledge, the direct numerical simulationsemployed in this work present the current state-of-the-artfor large scale simulations of electric interactions at thelevel of cell aggregates. We believe the simulation resultsprovide precise information about electric response of cellaggregates. In this manuscript we leverage the insightsdrawn from this direct numerical simulation data to de-velop and corroborate a theory for the time-dependentevolution of the TMP in cell aggregates under arbitraryapplied electric pulse.
D. Multiscale modeling strategy
Given the microscopic model of the electric interac-tions (6a)-(6e), we aim to infer effective theories at themulticellular level where tens of thousands of cells arepresent. In our work we focus on the effective propertiesof the aggregate using the effective medium theory , forexample see [67], accompanied with a dynamical modelfor the transient response of the system to an externalpulse using the
Fokker-Planck formalism [68, 69]. Be-low, we briefly describe and justify each component ofour modeling strategy. • Effective medium theory: Maxwell (1873) [37] wasthe first to study effective transport properties of a sta-tionary, random and homogeneous suspension of spher-ical particles dispersed in a background medium of uni-form conductivity. Maxwell’s assessment of effective con-ductivity was accurate to order O ( φ ) ( φ is the volumefraction of particles) and relied on his observation thatchanges in effective conductivity of a suspension of par-ticles was due to the average dipole moment of parti-cles. Exactly a century later, Jeffrey (1973) [70] ex-panded Maxwell’s estimation to order O ( φ ) using thegeneral method of Batchelor (1972) [71, 72] by consid-ering pairwise interactions between spherical particles.Batchelor’s work was focused on studying the effects ofhydrodynamic interactions between particles moving atlow velocity through a fluid on the effective viscosity of asuspension of particles. A remarkable result of his workwas finding the second order correction term to Einstein’sresult for the effective viscosity of a suspension of diluteparticles by systematically considering pairwise hydro-dynamic interactions between freely-moving spheres in alinear flow field.As part of the work presented here, we apply Batch-elor’s approach to the problem of conductive shelledspheres instead of the “homogenization” method thatemerged in the 80’s and is more commonly used in simi-lar application domains, e.g. in cardiac electrophysiology[73, 74]. Our choice is motivated by several reasons: (i)as was pointed out by Hinch (2010) [75], homogeniza-tion techniques are limited to periodic microstructures,which makes them less applicable for modeling finite sizemulticellular systems like tumor spheroids, (ii) present-ing this problem in Batchelor’s formalism allows for ap-plication of several existing results such as the influenceof the particles’ shape and arrangement on the effectiveconductivity of cell aggregates (even at maximum pack-ing fractions with touching spheres), the considerationof both near- and far-field interactions between particles,as well as the effects due to higher-order multipole mo-ment interactions of particle polarizations on the overallconductivity; we refer the interested reader to Batche-lor (1974) [76], Bonnecaze & Brady (1990) [77] and thereferences therein for more details. Lastly, (iii) Batch-elor’s method is based on ensemble averages of interac-tions among dispersed particles, that implicitly assumesindistinguishability of particles, which is in line with the Fokker-Planck formalism that we present next. • Fokker-Planck formalism: Cells in an aggregate areheterogeneous in shapes and electrical properties. There-fore, under an applied electric field, the evolutionary pathof a cell’s polarization is different from that of other cells,necessitating a probabilistic description of the inducedpolarizations. To this end, we describe the state of themulticellular system with the probability density of in-duced dipole moments on membranes and cytoplasms.We then derive a Fokker-Planck equation (FPE) thatdescribes the time evolution of the state-space probabil-ity density in response to an electric pulse, as well asthe many non-thermal small disturbances that influencethe states. Therefore, the Fokker-Planck equation notonly provides the stationary state of the system, but alsopredicts its dynamics far from equilibrium. The Fokker-Planck equation was first used by Fokker (1914) [68] andPlanck (1917) [69] and have been used to describe numer-ous systems such as the statistics of laser lights, or the ro-tations of dipole moments under various potentials. Thelatter was carried out by Debye who derived a Fokker-Planck equation for the rotational dynamics of polarizedmolecules and is used to describe the γ -dispersion dis-cussed in section I A. We refer to Risken (1984) [78] fora standard exposition of this topic, and to [79] for anoverview of applications in sciences and engineering.The basic idea of our treatment is to modify the bound-ary conditions on the cell membranes in order to add ap-propriate disturbances to the system parameters, whichin turn provides a Langevin equation for the induceddipole moment. Then, we transform the Langevin equa-tion to its corresponding Fokker-Planck partial differen-tial equation. Finally, we reduce the governing FPE byusing a moment-based approach and derive a simple setof ordinary differential equations (ODEs) that capturesthe evolution of the average and of the variance of in-duced polarizations. Importantly, the set of ODEs is sim-ple enough that it can be used for real-time predictionsand control of transmembrane potentials under arbitraryexternal electric stimulations and system parameters.The plan of this manuscript follows: in section II weuse Green’s theorem to decompose cellular polarizationinto its different components and compute effective con-ductivity of the medium. In section III we use thecontinuity of flux across cell membranes to develop aLangevin equation for the evolution of membrane polar-ization, thereafter we perturb the physical parameters inthis model and after standard averaging procedure weobtain the corresponding Fokker-Planck equation. Weprovide an analytic treatment of the governing FPE andleverage a moment based approach to compute the re-duced order ODE system for the statistical moments ofthe induced polarization density. Also, we argue in favorof a fractional order for time derivative in the FPE. Fi-nally, in section IV we present numerical results of ourmodel both in the time domain and frequency domain.We conclude our work in section V. II. COARSE-GRAINED REPRESENTATION
Our strategy is to use a multipole generator scheme torepresent the surface current density on individual cellmembranes in terms of an equivalent set of dipole mo-ments that reproduce an identical potential in a homog-enized medium.
A. Cells as arrays of layer potentials
We denote by u( r ), the electric potential at any point r within the aggregate; u( r ) satisfies Laplace’s equation.We treat this problem by viewing a membrane as a dipolelayer ( cf. see chapter 1 of [80] for more details) that isformed by two infinitesimally close surfaces with oppositecharge densities, as depicted in figure 1. In this work weonly consider a passive environment, i.e. that is freefrom current source or sink terms on cell membranes,and denote the passive current by J m . In this case, wecan relate the normal current density passing throughthe membrane by J ± m · n (where ± refers to the externaland internal side of the membrane) to the potential jumpacross membrane by writing − J ± m · n = σ c ∂ n u c = σ e ∂ n u e , (7)Similar to Geselowitz [81], we use Green’s theoremto treat this problem in terms of source current densi-ties and applied electric fields surrounding the domain.Consider two well-behaved functions ψ and φ definedinside and outside of cells and define the vector field F = σψ ∇ φ such that it is a continuous function of po-sition in enclosed volumes between boundaries. Substi-tuting F in Gauss’s theorem and wrapping thin mem-branes with discontinuity in between two close surfaces,one obtains Green’s theorem for non-homogeneous mix-tures with discontinuity across internal boundaries (seepage 53 of [82]), q (cid:88) j =1 (cid:90) A j [ σ c ( ψ c ∇ φ c − φ c ∇ ψ c ) − σ e ( ψ e ∇ φ e − φ e ∇ ψ e )] · d A j + p (cid:88) j =1 (cid:90) A j σ ( ψ ∇ φ − φ ∇ ψ ) · d A j = (cid:90) V [ ψ ∇ · ( σ ∇ φ ) − φ ∇ · ( σ ∇ ψ )] dv where q is the number of surfaces across which σ is discon-tinuous, p is the number of additional boundaries includ-ing macroscopic boundaries where electrodes are located,and V is the volume enclosed between Γ and all innersurfaces A j excluding surfaces of discontinuity. Here weadopt a convention that d A j is the measure of an areaelement of the surface A j that always points into theextra-cellular matrix.Geselowitz, in his case I , considers ψ = 1 /r and φ = u with r being the distance between surface or volume elements to any arbitrary point in the domain. Then itis straightforward to show that for any observation point x in the aggregate, we have4 πu ( x ) = (cid:88) m (cid:90) A m ( σ e u e − σ c u c ) ∇ (cid:48) ( 1 r ) · d A + (cid:90) Γ (cid:18) J b σ e r + u b ∇ (cid:48) r (cid:19) · d A , where Γ is the surface of the electrodes, u b is the electricpotential applied at the electrodes, J b = − σ e ∇ u , and A m is the surface of membranes that points into the extra-cellular matrix. Moreover, as in case II of Geselowitz, wecould let ψ = 1 /r and σφ = u to show that for infinites-imally thin membranes the electric potential is given by:4 πu ( x ) = (cid:88) m (cid:90) A m (cid:20)(cid:18) σ e − σ c (cid:19) J ± m r + [ u ] ∇ (cid:48) ( 1 r ) (cid:21) · d A + (cid:90) Γ (cid:18) J b σ e r + u b ∇ (cid:48) r (cid:19) · d A , (8)Hence, using the divergence theorem one can obtain thealternative formulation:4 πu ( x ) = (cid:88) m (cid:90) A m [ u ] ∇ (cid:48) ( 1 r ) · d A m + (cid:88) c (cid:90) V c (cid:18) σ e − σ c (cid:19) J · ∇ (cid:48) ( 1 r ) dV + (cid:90) Γ (cid:18) J b σ e r + u b ∇ (cid:48) r (cid:19) · d A , (9)where r = | x − x (cid:48) | , ∇ (cid:48) = ∂/∂ x (cid:48) and the summation isover the membranes of the enclosed cells within Γ. Thefirst term on the right-hand side describes the influenceof the polarized membranes, the second term capturesthe contribution from the cytoplasms, and the last termrepresents the influence of electrodes. We also empha-size that if the electrodes are not in direct contact withthe extra-cellular matrix (e.g. there is a gap with lowconductivity between Γ and the outer surface of the ag-gregate), one has to also include the extra contributionfrom the outer surface ( A o ), u o ( x ) = 14 π (cid:90) A o E o r · d A . In addition to the contribution from the electrodes, equa-tion (8) decomposes the electric potential at any point inthe volume as a superposition of a monopole/single layerand a dipole layer on each membrane. Remarkably, themembrane integral over transmembrane jump is analo-gous to the contribution from a surface current dipoledensity: u p ( x ) = (cid:88) m πσ e (cid:90) A m D ( x (cid:48) ) ∇ (cid:48) (cid:18) r (cid:19) · d A m where the current dipole surface density ( i.e. currenttimes distance between sources) is defined by: D ( x ) = σ e [ u ]Then, a point dipole on the membrane is expressed as δ P = σ e [ u ] d A m and one may model the induced trans-membrane potential as a resultant current dipole on eachcell.Thus, we approximate each membrane with surface A i by a resultant dipole of strength P i = (cid:90) A i σ e [ u ] d A , (10)and we call P the dipolar polarization .Furthermore, we observe that the membrane integralover the transmembrane current density resembles thecontribution from a monopole current layer , u s ( x ) = (cid:88) m π (cid:90) A m (cid:18) σ e − σ c (cid:19) J ± m r · d A m = (cid:88) m πσ e (cid:90) A m ( σ c − σ e ) E c r · d A m where E c refers to the electric field at the inner surface ofthe membrane. Here, we identify an induced polarizationdensity of ( σ c − σ e ) E c · n / (4 πσ e ) over cell membranes.Alternatively, this term in its volume integral form reads u s ( x ) = (cid:88) c πσ e (cid:90) V c ( σ c − σ e ) E · ∇ (cid:48) (cid:0) r (cid:1) dV, which can be interpreted as an ‘extra flux density’, de-noted by τ and is zero everywhere in the extra-cellularmatrix while in the cells is given by τ ( x ) = ( σ c − σ e ) E = J ( x ) + σ e ∇ u ( x ) . Hence u s ( x ) = (cid:88) c πσ e (cid:90) V c τ ( x (cid:48) ) · ∇ (cid:48) (cid:0) r (cid:1) dV, which again resembles the electric potential of a volumedipole density τ . Therefore, we define the instantaneouspolarization ( S ) to model the polarization of cell cyto-plasms: S i = (cid:90) V i τ dV, or, equivalently, in terms of the potential in the cytoplas-mic side of the membrane as S i = ( σ e − σ c ) (cid:90) A i u c d A . Then the net polarization can be defined by M i M i = P i + S i = (cid:90) A i ( σ e u e − σ c u c ) d A We also note that this result could be directly inferredusing Green’s theorem by letting φ = u and ψ = 1 /r , e.g. see equation 29 of Geselowitz 1967. Furthermore one canrelate the values of P and S through the equation: P σ e = (cid:90) A u e d A + S σ c − σ e (11)We observe that in general the dipolar and instantaneouspolarizations may have different orientations dependingon the symmetries of the exterior potential. To assessthis relation, we use Gauss’ theorem for a closed surfaceΓ enclosing N internal closed surfaces ( cf. see chapter IIIof Smythe, note the minus sign is due to our conventionthat the normal direction points into the extra-cellularmatrix), N (cid:88) j =1 (cid:90) A j u d A + (cid:90) Γ u d A = − (cid:90) V e ∇ u dV, (12)where V e is the volume of the extra-cellular matrix ex-cluding cells. Therefore, for the cell aggregate n < P >σ e = (cid:80) Nj =1 (cid:82) A j u d A V + n < S >σ c − σ e , (13)where, n < P > = (cid:80) Nj =1 P j V and n < S > = (cid:80) Nj =1 S j V , with n the number density of cells in the mixture. Notethat the volume fraction φ is related to the number den-sity n via φ = nV c . Using Gauss law (12) with equation13 we obtain n < P >σ e = n < S >σ c − σ e − V (cid:90) V e ∇ u dV − V (cid:90) Γ ud A . Furthermore, we define the applied electric pulse on theboundary, E ext , as: E ext ≡ V (cid:90) Γ ud A . We also recall that S is related to the volume averagedelectric field in cell cytoplasms by: φ ¯ E c ≡ − V cells (cid:88) j (cid:90) V j ∇ u dV = n < S >σ c − σ e , while the volume averaged external field is simply(1 − φ ) ¯ E e ≡ − V (cid:90) V e ∇ u dV = E ext + n < P >σ e + n < S >σ e − σ c and the volume average electric field inside the mem-branes is related to the dipolar polarization by1 V cells (cid:88) j (cid:90) V (cid:48) j ∇ udV = 1 V cells (cid:88) j (cid:90) Γ j [ u ] d A = n < P >σ e . Because the volume V is partitioned into three parts V = V e ∪ V c ∪ V (cid:48) c , where V (cid:48) c is the volume occupied by cellmembranes, we conclude that − < ∇ u > = E ext = φ ¯ E c + (1 − φ ) ¯ E e − n < P >σ e . In particular, we note that φV = (cid:80) j V j , and we definethe dipolar polarization per cell volume as p j = P j V j toobtain the effective dipole moment per cell volume,¯ p = (cid:80) j V j p j (cid:80) j V j , and we can write n < P > = φ ¯ p . B. Frequency domain model for cell dipoles
We consider the analytical solution of a spherical cellof radius R centered within a spherical domain of radius R under a Dirichlet potential E ( t ) R cos θ at the outerboundary (a local electric field E = − E k is consideredat the surface of a sphere of radius r = R from centerof the cell). In this case the membrane voltage satisfies C m ∂ [ u ] ∂t + ( S L − B )[ u ] = AE ( t ) R cos θ, where A = 3 σ c σ e R K and B = − σ c σ e ( R + 2 R R ) K, with K − = R ( σ e − σ c )+ R (2 σ e + σ c ) so the coefficientsare: A = 3 σ c σ e /R σ e + σ c + φ ( σ e − σ c ) , and B = − (2 + φ ) σ c σ e /R σ e + σ c + φ ( σ e − σ c ) , where the volume fraction of cells is given by φ = R /R .Furthermore, we define three independent parametersthat characterize the solution:˜ σ = 2 σ e + σ c + φ ( σ e − σ c ) ,η = 1 + S L R ˜ σ (2 + φ ) σ c σ e ,τ = ˜ σR C m (2 + φ ) σ c σ e . In the frequency domain the solutions are given by,[˜ u ] = 3 R φ · η + jωτ · ˜ E ( ω ) · cos θ ˜ u c = ˜ α c ˜ E ( ω ) · r · cos θ ˜ u e = (˜ α e r + ˜ β e r ) · ˜ E ( ω ) · cos θ with ˜ α c = 3 σ e ˜ σ · (cid:18) η − jωτη + jωτ (cid:19) ˜ α e = (cid:18) σ c + 2 σ e ˜ σ − σ c ˜ σ · φ φ · η + jωτ (cid:19) ˜ β e = (cid:18) σ e − σ c ˜ σ + 3 σ c ˜ σ ·
12 + φ · η + jωτ (cid:19) · R The external and internal electric fields are computedgiven the electric potential ( r = r ˆ r )˜ E e ( ω ) = ˜ α e ˜ E ( ω ) −
3( ˜ β e ˜ E ( ω ) · ˆ r )ˆ r − ˜ β e ˜ E ( ω ) r , (14)and˜ E c ( ω ) = ˜ α c ˜ E ( ω ) , (15)which corresponds to a uniform external field superim-posed by the electric field of a net dipole moment˜ M ( ω ) = − πσ e ˜ β e ˜ E ( ω ) . (16)It can be easily verified that indeed u e ( R ) = ER cos θ as expected. We emphasize that we only impose the tan-gential component of electric field at r = R and not itsradial component, as evident by equation (14).As we discussed before, we represent this solutionby defining dipole moments over membranes and cyto-plasms. It is straightforward to compute the dipole mo-ments from the basic definitions of the previous section,which lead to˜ P ( ω ) V c = − σ e φ · η + jωτ · ˜ E ( ω ) (17)and˜ S ( ω ) V c = − σ e ( σ e − σ c )˜ σ · η − jωτη + jωτ · ˜ E ( ω ) , (18)where V c is the volume of a cell. One can verify thatindeed we have ˜ P + ˜ S ≡ ˜ M from equation (16). More-over, it is straightforward to verify that ˜ P and ˜ S arerelated through equation (11). We note the minus signin dipole moments stems from our mathematical defini-tion for jump in solution across membrane, that is thevalue of solution in the exterior minus that of the cyto-plasm, which is opposite the usual convention in biology.Therefore, in making the connection with other works,care must be taken in using consistent signs at this step.Note also that the influence of the membrane conductiv-ity is captured in the parameter η . For example, for aninsulating membrane where σ m (cid:28) σ c , we have η ≈ i.e. S = 0.Equations (17)–(18) enable the definition of the cellu-lar polarizability coefficients P = σ e α p E , S = σ e α s E ,and M = σ e α E with α p = −
32 + φ · η + jωτ ,α s = − σ e − σ c )˜ σ · η − jωτη + jωτ , and α = α p + α s . Also note that ˜ α e is related to the electric polarizability α p via ˜ α e = 2 σ e + (1 + φα p ) σ c ˜ σ . The average polarization of the whole aggregate is givenby: n < ˜ P > = φ · α p · σ e ˜ E (19) n < ˜ S > = φ · α s · σ e ˜ E (20)Importantly, ˜ E can be related to the applied pulse ˜ E ext by averaging equation (14) within a spherical shell vol-ume between the membrane and R . By integration,the dipolar contribution is zero, and we establish a re-lationship with the average external electric field, ¯ E e ,˜ α e ˜ E = (1 − φ ) ¯ E e . We already showed that(1 − φ ) ¯ E e = ˜ E ext + n < ˜ P >σ e + n < ˜ S >σ e − σ c , therefore the local electric field is given by˜ E = κ − ˜ E ext , (21)with κ = α e − φα p − φ σ e α s σ e − σ c , (22)which simplifies to κ = (2 + 3 φ ) σ e + σ c ˜ σ − φ σ c (2 + φ )˜ σ ( η + jωτ ) . (23)For example figure 3 illustrates the magnitude of E fortwo different membrane conductivities.At the aggregate level, we can define the membranesusceptibility, n < P > = σ e χ p E ext , the cytoplasm sus-ceptibility, n < S > = σ e χ s E ext , as well as the overall cellsusceptibility, n < M > = σ e χ E ext , that relate the ap-plied electric pulse to the induced polarization densities.Equation (21) allows us to relate the applied pulse to theinduced polarizations, therefore we obtain the suscepti-bility coefficients χ p = φα p κ and χ s = φα s κ , which also imply that(1 − φ ) ¯ E e E ext = α e κ and φ ¯ E c E ext = σ e χ s σ c − σ e . The time-domain version of E can be found by express-ing the complex factor η + jωτ in equation (23) in termsof < P > . By noting that η + jωτ = − φσ e ˜ E (2 + φ ) n < P > , we obtain that κ = (2 + 3 φ ) σ e + σ c ˜ σ + φ σ c ˜ σ n < ˜ P >σ e ˜ E , which, substituted into equation (21) leads to E ( t ) = ˜ σ E ext ( t ) − νφn < P ( t ) > (2 + 3 φ ) σ e + σ c where ν = σ c /σ e . For the purpose of developing time-domain equations it is also useful to write the cytoplasmpolarization in terms of the membrane polarization: n < S > = − φ σ e − σ c ˜ σ (cid:18) σ e E + 2 + φ φ n < P > (cid:19) . C. Effective conductivity
Maxwell (1873) [37] first considered the problem ofcalculating effective conductivity coefficients for dilutespherical inclusions. A century later, Jeffrey (1973) [70]included pairwise interactions into Maxwell’s theory, forincreasing the validity of the estimated effective con-ductivity for higher concentrations. Chiew and Glandt(1983) [83] considered a more realistic pair-correlationfunction to improve on the accuracy of Jeffrey’s result.In parallel, Hasselman and Johnson 1987 [84] includedthe effect of interfacial resistance to Maxwell’s theory,which was subsequently integrated with Jeffrey’s theoryby Chiew and Glandt [85]. In this section, we derive theeffective conductivity based on the work of Batchelor fortransport phenomena in two-phase media composed ofan statistically homogeneous suspension of particles withrandom configurations (see e.g. [76, 77, 86]). We empha-size that for improved accuracy one has to modify thisapproach to include divergences of all higher order multi-pole moments; we refer the interested reader to [87, 88].First, we define the average flux in a volume largeenough to include many cells, < J > = 1 V (cid:90) V J dV, and we seek a linear relationship between the average fluxand the potential gradient < J > = − ¯ σ < ∇ u >, (24)0 FIG. 3. The magnitude of the local electric field at different frequencies and volume fractions for ( σ e , σ c ) = (1 . , .
6) [
S/m ]and S L = 1 . S/m ] (left) and S L = 1 . × [ S/m ] (right). where the proportionality coefficient defines the effectiveconductivity . We decompose < J > into three differentregions, < J > = 1 V (cid:90) V − (cid:80) i V i J dV + 1 V (cid:88) i (cid:90) V i ∪ V (cid:48) i J dV = 1 V (cid:90) V − σ e ∇ u dV + 1 V (cid:88) i (cid:90) V i ∪ V (cid:48) i τ dV where (cid:80) i V i is the volume occupied by cells. Last expres-sion is obtained by replacing J k = − σ k ∇ u ≡ − σ e ∇ u + τ k such that τ k = ( σ e − σ k ) ∇ u . In the membrane we ap-proximate ∇ u = [ u ] n /h , therefore, < J > = − σ e < ∇ u > + n < S > + (cid:0) − σ m σ e (cid:1) n < P > (25)which in terms of the dipolar polarization and the exter-nal electric field is given by < J > = σ e E ext (cid:18) ν + 3 φν ν + 3 φ (cid:19) + n < P > (cid:18) − σ m σ e − σ e − σ c ˜ σ [2 + φ − νφ φ + ν ] (cid:19) . Moreover, using the definition of the effective conductiv-ity (24) and the fact that − < ∇ u > = E ext , we deducethat the parallel ( (cid:107) ) and the transverse ( ⊥ ) componentsof the effective conductivity are given by¯ σ (cid:107) σ e = < J (cid:107) >σ e E ext and ¯ σ ⊥ σ e = < J ⊥ >σ e E ext . Particularly, in the frequency domain, the parallel com-ponent satisfies¯ σσ e = 1 + χ s + (1 − σ m σ e ) χ p . (26)For infinitely conductive membranes, where η → ∞ ,equation (26) can be simplified by noting that σ m σ e = hR · (2 + φ ) σ c ˜ σ · ( η −
1) (27)and that the membrane polarization vanishes. In thiscase, equation (26) reduces to Maxwell’s equation for theeffective conductivity of a dilute suspension [89]:¯ σσ e → − φ − ν + νh/R ν + 3 φ (28)and h/R →
0. Note that Maxwell’s result is only a O ( φ )estimate of the effective conductivity, and that, to itslimit of validity, equation (28) coincides with Maxwell’sapproximation (see e.g. [70]).On the other hand, one could identify relative complexpermittivity (cid:15) ∗ = (cid:15) (cid:48) − j(cid:15) (cid:48)(cid:48) through its definition, i.e. thecurrent J is related to an alternating applied field E via J = σ e E + jω(cid:15) (cid:15) ∗ E . Comparing with equation (25), wecan write jω(cid:15) (cid:15) ∗ E ext = (1 − σ m σ e ) n < P > + n < S >, with (cid:15) ∗ = σ e jω(cid:15) · ((1 − σ m σ e ) χ p + χ s ) = (cid:15) (cid:48) − j(cid:15) (cid:48)(cid:48) , (a) (b)(c) (d) FIG. 4. Equation (30) for ν = 0 .
46 with σ e = 1 . S/m ], σ c = 0 . S/m ], and the cell radius set to 7 [ µm ]. (a,b) Low membraneconductance S L = 1 . S/m ], (c,d) high membrane conductance S L = 1 . × [ S/m ]. where (cid:15) (cid:48) is the dielectric constant, (cid:15) (cid:48)(cid:48) is the loss factor ofthe material and (cid:15) = 8 . × − [ F/m ] is the dielec-tric permittivity of free space. Furthermore, the complexadmittance is given by Y ∗ = ( σ e + jω(cid:15) (cid:15) ∗ ) A el /H . Theadmittance being the inverse of the impedance, Z , wehave: Z = HA el · σ e + jω(cid:15) (cid:15) ∗ = HA el · E ext · k ( σ e E ext + (1 − σ m σ e ) n < P > + n < S > ) · k (29)where H is the distance between the electrodes and A el is the surface area of one electrode. We designate Z e = H/ ( σ e A el ) to define the dimensionless impedance as ZZ e = σ e E ext · k ( σ e E ext + (1 − σ m σ e ) n < P > + n < S > ) · k . = 11 + (1 − σ m σ e ) χ p + χ s Lastly, the complex conductivity σ ∗ of the material isrelated to the admittance according to Y ∗ = σ ∗ A el /H =1 /Z which implies that: σ ∗ = σ e · (1 + (1 − σ m σ e ) χ p + χ s ) . Therefore, the values for admittance and the effectiveconductivity coincide.2
Effects of cellular pairwise interactions.
Cell-cell in-teractions influence the effective conductivity of the ag-gregate, for example in the case of spherical particlesJeffrey (1973) [70] showed that pairwise interactions pro-duce a correction term of order O ( φ ) to the effectiveconductivity of a random dispersion of spherical parti-cles. Importantly, Jeffrey used the twin spherical har-monics invented by Ross (1968) [90] to account for two-particle interactions. This method was later applied tocoated spheres by Lu and Song (1996) [91]. [91] derivedthe general expression for the effective conductivity of arandom suspension of coated spheres, that is accurate upto order O ( φ ):¯ σσ e = 1 + 3 φθ + 3 φ θ − φθ + K ∗ φ − φθ , (30)where 3 φθ ≡ χ s + (1 − σ m σ e ) χ p is the polarizability fac-tor. K ∗ accounts for higher order interactions due to de-tailed pair distribution of particles; here we neglect thislast term as we are not considering detailed informationabout the microstructure of the aggregate (see Hasselmanand Johnson (1987) [84] for a similar result). As shown by[91], this estimate stays within Hashin-Shtrikman bounds[92] up to high volume fractions of about φ ≈ .
6. Figure4 illustrates the dependence of the effective conductivityon the volume fraction and on the frequency predictedby equation (30) after setting K ∗ = 0. III. COARSE-GRAINED DYNAMICS
In the present modeling approach, cell-level dipole mo-ments are the resolved observables that relate cellularproperties to multicellular features. In this section, weintroduce time-domain governing equations for the dipolemoments.We consider a spherical cell of radius R immersed in amean electric flux < J > , C m ∂ [ u ] ∂t + (cid:18) S m + (2 + φ ) σ e σ c R ˜ σ (cid:19) [ u ] = 3 σ c σ e ˜ σ E cos θ, (31)and we seek a governing equation for the membranedipole by integrating equation (31) over cell membranesand multiplying by σ e , C m ddt P + (cid:18) S m + (2 + φ ) σ e σ c R ˜ σ (cid:19) P = − σ e σ c ˜ σ A E . Here we assumed a uniform conductance over the cellmembranes. We divide both sides with the cell volumeand obtain: C m ddt p + (cid:18) S m + (2 + φ ) σ e σ c R ˜ σ (cid:19) p = − σ e σ c R ˜ σ E Due to this mean-field approximation, each cell evolvesostensibly independently from each other. Therefore, we define the coarse grained electrodynamics of an individ-ual cell with ˙ p = − γ p − α u ( t ) , (32)where we represent the time-dependent model for theelectric flux with u ( t ) and define α = 3 σ e σ c C m R ˜ σ γ = S m C m + σ e σ c RC m ˜ σ (2 + φ )and the stimulating field is that of the mean electric fieldin the matrix at any given time t , u ( t ) ≡ σ e E ( t ) = σ e κ − E ext ( t ) . Note that u ( t ) is the mean field approximation for theexternal electric current that each cell feels.Starting from the microscopic equation (32) we shallderive a mesoscopic model for an ensemble of cells. Thebasic idea is to view an ensemble of dipoles as ran-dom variables, and subsequently treat equation (32) asa Langevin equation for the dynamical evolution of therandom variables. In 1908, Langevin introduced the con-cept of equation of motion of a random variable [93] andthrough his formulation of the dynamical theory of Brow-nian motion, he initiated the subject of stochastic differ-ential equations [94].To bridge the microscale to the mesoscale, we are in-terested to know the probability distribution of dipolemoments in an aggregate. In stochastic systems ( e.g. in many condensed matter systems that are in contactwith a heat bath) a successful strategy is to start froma Langevin equation describing the evolution of a sin-gle particle. Then, through appropriate averaging proce-dures, one arrives at a Fokker-Planck equation describingthe evolution of the probability distribution of that par-ticle. Eventually, the independence assumption impliedin the mean field approximation allows to define the to-tal probability distribution W ( { p k } , t ) as the product ofthat of individual particles W ( { p k } , t ) = Π i W i ( p i , t ).Unfortunately, this procedure fails in the system of cellaggregates due to the deterministic nature of the elec-trodynamics of cells. Through direct numerical simula-tions we know that the evolutionary trajectory of celldipoles is not a stochastic process; in fact equation (32)already suggests that polarizations are given by a deter-ministic equation. In the next subsection, we solve thisproblem by considering the randomness in cell param-eters (see figure 5) and devise a stochastic interpreta-tion. We also mention the interesting work of Takayasu et al. (1997) [95] who took a similar strategy to analyzethe conditions for the emergence of power-law distribu-tions in specific discrete Langevin models of the form x ( t + 1) = b ( t ) x ( t ) + f ( t ) where both b ( t ) and f ( t ), arerandom variables.3 FIG. 5. Distribution of α and γ parameters used in numericalsimulation. The densities follow exponential profile in themiddle range. A. The indistinguishable stochastic replica
The main source of randomness in the dielectric re-sponse of multicellular systems is the diversity of cell pa-rameters, i.e. namely α and γ . Here we exercise analternative viewpoint to replace a diverse ensemble ofdeterministic cells with an indistinguishable ensemble ofstochastic elements. In particular we consider a fiducial random-walk process in cell parameters such that the em-pirical probability density of cells matches that of actualones at any timestep. This is achieved by defining tworandom processes A ( t ) and B ( t ) for each cell such that α i = ¯ α + A i ( t ) ,γ i = ¯ γ + B i ( t ) , and ¯ γ = < γ >, ¯ α = < α > . Figure 5 illustrates the distribution of these parameters.It is evident that around their mean values the distribu-tions follow an exponential profile (note that the centralpart of this figure is linear and the y-axis is in logarith-mic scale), however as a first step we approximate these random processes with Gaussian white noise < B i ( t ) > = 0 < A i ( t ) > = 0 < A i ( t ) A j ( t (cid:48) ) > = α (cid:48) δ ij δ ( t − t (cid:48) ) < B i ( t ) B j ( t (cid:48) ) > = γ (cid:48) δ ij δ ( t − t (cid:48) ) < A i ( t ) B j ( t (cid:48) ) > = (cid:15)α (cid:48) γ (cid:48) δ ij δ ( t − t (cid:48) )where | (cid:15) | ≤ p i = − ¯ γ p i − ¯ α u ( t ) − p i B i ( t ) − u ( t ) A i ( t ) . (33)Equation (33) is a Langevin model with both additiveand multiplicative noise terms. In analogy, the additivenoise term models a heat bath acting on the dipole and themultiplicative noise term models effects of a fluctuatingbarrier. Interestingly, such hybrid models of arithmeticand geometric Brownian motions have many applicationsin Physics [96, 97], Biology [98] and finance [99].We take equation (33) along each spatial dimensionas an independent stochastic differential equation in therandom variable X t ≡ p k , then we combine the two Brow-nian motions into a single Brownian motion to arrive atthe stochastic differential equation dX t =( − ¯ αu − ¯ γX t ) dt + (cid:113) α (cid:48) u + 2 (cid:15)γ (cid:48) α (cid:48) uX t + γ (cid:48) X t dW t (34)where W t is a Wiener process. This is in fact a memberof Pearson diffusions that were extensively considered byForman and Sørenson [100].Here we emphasize that the replica stochastic aggre-gate must render a similar distribution of polarizationsas in the actual system, therefore in general the processes A ( t ) and B ( t ) are not Markov processes with Delta au-tocorrelation in time. In fact the specificities of theserandom processes must be tuned to match the actual dis-tribution. In the next section we simply consider whitenoise as a first approximation in this direction. B. The Fokker-Planck equation
Equation (34) constitutes the Stratonovich vectorstochastic differential equations of a Langevin equationwith a multiplicative noise term. The generic form ofthe coupled stochastic differential equation in terms ofstochastic variables ξ i ’s reads d { ξ } i dt = h i ( ξ , ξ , ξ , t ) + g ij ( ξ , ξ , ξ , t ) λ j ( t )with g ij λ j are multiplicative noise terms that producethe noise-induced drift and diffusion components. λ hasthe following property, < λ i ( t ) > = 0 , < λ i ( t ) λ j ( t + τ ) > = δ ij δ ( τ )4Then these equations can be treated as the starting pointfor deriving the corresponding Fokker-Planck equation.To this end, let ξ = p x , ξ = p y , and ξ = p z and W ( ξ , ξ , ξ , t ) dξ dξ dξ be the probability of finding adipole in dξ dξ dξ at time t , then the Fokker-Planckequation reads ∂W∂t = − ∂∂ξ i ( D i W ) + 12 ∂ ∂ξ i ∂ξ j ( D ij W )where the Einstein’s summation rule is implied. By useof the Langevin equations, one can evaluate the statis-tical averages < · > in order to get the following set ofequations for the drift and the diffusion coefficients [101]: D i = h i ( { ξ } , t ) + 12 g kj ( { ξ } , t ) ∂∂ξ k g ij ( { ξ } , t ) ,D ij = g ik ( { ξ } , t ) g jk ( { ξ } , t ) . Hasegawa (2008) [102, 103] considered solutions to theFokker-Planck equations associated with Langevin equa-tion (33) in the Stratonovich stochastic calculus; further-more, Mortensen (1979) [104] derived the Fokker-Planckequation associated to SDE 34 according to Ito stochas-tic calculus. Fortunately this system is separable andwe can treat it one dimension at a time for our analysis, i.e. the other two dimensions can be identically treated.The governing FPE for the variable x ≡ p k with theStratonovich interpretation reads: ∂∂t W ( x, t ) = ∂∂x (cid:20) ¯ γx + ¯ αu ( t ) − χ γ (cid:48) x + (cid:15)γ (cid:48) α (cid:48) u ( t )) (cid:21) W ( x, t )+ 12 ∂ ∂x (cid:20) γ (cid:48) x + 2 (cid:15)γ (cid:48) α (cid:48) u ( t ) x + α (cid:48) u ( t ) (cid:21) W ( x, t ) (35)where χ = 0 , C. Analytical Treatment
The invariant probability distribution has a densitythat satisfies: ddx W ( x ) = − (¯ γ + χ γ (cid:48) ) x + (¯ α + χ (cid:15)α (cid:48) γ (cid:48) ) u γ (cid:48) x + (cid:15)α (cid:48) γ (cid:48) ux + α (cid:48) u W ( x ) (36)Equation (36) resembles the invariant density of Pearsondiffusion processes that are characterized with a lineardrift and quadratic diffusion coefficients. Therefore, wefind that the Fokker-Planck equation (35) falls in the cat-egory of Pearson diffusion processes [105], whose station-ary probability density is invariant under translation andscale transformations. Statistical properties of this classof models have been analyzed by Forman and Sørenson[100] (for a brief summary see section 1.13.12 of [106]). Fundamentally, Pearson diffusions are viewed as the so-lutions to the following stochastic differential equation inthe canonical parameterization: dX t = − θ ( X t − ˆ µ ) dt + (cid:113) θ ( a X t + b X t + c ) dW t , with θ > a , b , c are shape parameterssuch that the diffusion coefficient is well defined, and W t is a Wiener process.Pearson processes can lead to a variety of distributionsdepending on the parameters of the drift and the diffusioncoefficients such as heavy or light tailed and symmetric orskewed profiles [100]. We briefly report six basic subfam-ilies from [100, 106] that are determined using criteria onthe degree of the diffusion polynomial in the denomina-tor of equation (36) denoted here by deg , the sign of theleading coefficient (that in our case is strictly positive),and the discriminant ∆ = b − ac :1. if deg = 0 : A Ornstein-Uhlenbeck process with anormal invariant density.2. if deg = 1 : If 0 < ˆ µ ≤ µ >
1, we obtain a gammainvariant density.3. if deg = 2 and ∆ > and a < : A Jacobi diffusionwith a Beta invariant density.4. if deg = 2 and ∆ > and a > : A Fisher-Snedecor process with a Fisher-Snedecor invariantdensity.5. if deg = 2 and ∆ = 0 : A Reciprocal gamma pro-cess with an inverse gamma invariant density.6. if deg = 2 and ∆ < : for ˆ µ (cid:54) = 0 we obtain aStudent diffusion with a skewed t invariant density,while for ˆ µ = 0 we obtain a scaled t -distribution.In the current case, and under a non-zero applied pulse,we can establish that the discriminant is always negative∆ = − (cid:0) α (cid:48) γ (cid:48) u (cid:1) (1 − (cid:15) ) < . Therefore, we expect that the probability density of thedipole moments along the z-axis (parallel to the appliedpulse) is best described by a skewed Student distribu-tion (also known as Pearson type IV distribution). Eventhough in the transverse direction the mean electric pulseis negligible based on our mean field model, we can not to-tally neglect the influence of the electric fluctuations. Tofirst order approximation we treat the transverse direc-tion by setting ˆ µ = 0 while preserving the same diffusionterm as in the parallel direction. This also ensures a sym-metric probability density for positive or negative values.Therefore, we conclude that the distribution in the trans-verse direction must follow a scaled t -distribution (alsoknown as Pearson type VII distribution).5
1. Stationary Probability Density a. In the direction parallel to the applied pulse.
Wecan directly integrate equation (36) to identify the invari-ant distribution W s ( x ; ν, c, a, λ ) = K exp (cid:8) c tan − (cid:0) x − λa (cid:1)(cid:9)(cid:0) (cid:0) x − λa (cid:1) (cid:1) ν , (37)where a = α (cid:48) uγ (cid:48) (cid:112) − (cid:15) , λ = − (cid:15)α (cid:48) uγ (cid:48) ,ν = χ γγ (cid:48) , c = (cid:15)α (cid:48) ¯ γ − ¯ αγ (cid:48) α (cid:48) γ (cid:48) √ − (cid:15) , and K = (cid:12)(cid:12) Γ( ν + ic ) (cid:12)(cid:12) a √ π Γ( ν )Γ( ν − / . Equation (37) provides the average and the variance ofthe polarizations in the stationary state ( cf. original workof [107], and [108] for a useful guide): µ ≡ E ( X t ) = acν − λ = (cid:15)α (cid:48) γ (cid:48) u − ¯ αu ¯ γ − γ (cid:48) , (38) σ ≡ E ( X t ) − E ( X t ) = a [( ν − + c ]( ν − (2 ν − γ (cid:48) µ + 2 (cid:15)γ (cid:48) α (cid:48) uµ + α (cid:48) u γ − γ (cid:48) ) , (39)where we have let χ = 1 to obtain the last equalities.Equation (37) can be independently verified by compar-ing it to model B of Hasegawa [103] via the change ofnotations λ H ≡ ¯ γ , I H ≡ − ¯ αu , α H ≡ γ (cid:48) , c H ≡ c , b H ≡ ν , f H ≡ − λ and β H ≡ α (cid:48) u (the H subscript indicatesHasegawa’s notation). Moreover, Hasegawa [103] derivedapproximate equations for the evolution of the averageand the variance in this model that we report here forcompleteness: dµdt = − ¯ γµ ( t ) − ¯ αu + γ (cid:48) µ ( t ) + (cid:15)γ (cid:48) α (cid:48) u, (40) dσ dt = − γ − γ (cid:48) ) σ ( t ) + γ (cid:48) µ ( t ) + 2 (cid:15)γ (cid:48) α (cid:48) uµ ( t ) + α (cid:48) u , (41)which provide the analytical formula for the evolution ofthe probability density throughout the polarization pro-cess when replacing the static values ν and c in equation(37) by their dynamic counterparts: ν ( t ) = a + ( µ ( t ) − λ ) + 3 σ ( t ) σ ( t ) ,c ( t ) = (cid:18) a + ( µ ( t ) − λ ) + σ ( t ) aσ ( t ) (cid:19) ( µ ( t ) − λ ) . Figure 9 gives a comparison between the results obtainedwith model (41) and the direct numerical simulation of Mistani et al. [59]. A few remarks follow: (i) Hasegawa’smoment equations are valid approximations up to or-der O (( δx ) ) about the mean value of dipolar polar-ization per cell volume, and (ii) in the current model,the dynamics for µ ( t ) appears to be decoupled from σ ,which is merely the result of assuming a linear depen-dence for the conductance term F ( x ) ≡ − ¯ γx , as wellas a linear assumption for the multiplicative noise factor G ( x ) ≡ x . Breaking either of these assumptions intro-duces extra contributions from the variance to the meandynamics. However, such modifications would alter thestationary distribution of polarizations, which the cur-rent model captures well. Therefore we expect that theobserved discrepancies in the temporal evolution betweenthe model prediction and the direct numerical simulationmost likely stem from the nature of the noise term. Thiscould be alleviated by introducing fractional-order tem-poral derivatives. Indeed, figure 9 actually proves thatthe transient dynamics of dipolar moments does not sim-ply follow an exponential function. b. In the direction perpendicular to the applied pulse. In the transverse direction, the stationary density fol-lows a symmetric t -distribution, which is obtained bysetting the skewness parameter to zero in the skewed t -distribution (37), i.e. set c = 0 and all previous equationsare valid. This indicates that the following identity reg-ulates the transverse diffusion, (cid:15) ⊥ α (cid:48)⊥ ¯ α = γ (cid:48)⊥ ¯ γ
2. Statistical moments
In the type IV Pearson diffusions considered in thiswork, for n ≥ n th statistical moment can be com-puted using the recurrence formula [108]: µ n = a ( n − ν − [2( ν − − ( n − × (cid:26) c ( ν − µ n − + a (( ν − + c ) µ n − (cid:27) (42)where by definition µ = 1 and µ = 0. It is straight-forward to check the consistency between equations (42)and (39).An important notion is the regime of existence for eachof the moments. The first moment exists for ν > < x > = ±∞ . The variance exists for ν > / ν > ν > /
3. Fitting statistical moments
We use the simple moment fitting approach intro-duced by Karl Pearson [107] (also see Heinrich’s ex-cellent guide [108]) to infer the four model parameters6
FIG. 6. (Left) The maximum values of the statistical moments diverge for larger k and (right) the evolution of the second tothe sixth statistical moments of p z [A / mm ] in direct numerical simulations, i.e. m k = < ( p z − < p z > ) k > .FIG. 7. Evolution of dipole moments using direct numericalsimulation. ( ν, c, a, λ ), which characterize the stationary probabilitydensity (37). Given simulation data, we can directly mea-sure the first four statistical moments using the recur-rence relation in equation (42), i.e. we directly calculate < x > , µ , µ and µ . Afterwards, we can infer theunknown parameters using mean, variance and some in-termediate quantities defined via the third and fourthmoments: (cid:112) β ≡ µ µ / = 2 cν − (cid:115) ν − ν − + c β ≡ µ µ = 3(2 ν − ν + 2)(( ν − + c ) − ν − ]( ν − ν − ν = 5 β − β − β − β − ,c = ( ν − ν − √ β (cid:112) ν − − β ( ν − ,a = (cid:114) µ [(2 ν − − β ν − ] ,λ = < x > − √ µ β ( ν − . Figure 8 illustrates the fitted model and its parameters.We observe that the model given in (37) perfectly de-scribes the results provided by our direct numerical sim-ulation. We compare the predictions of our model withthe dynamics of the average and of the variance of dipolemoments from our direct numerical simulations in fig-ures 9–10. A few observations regarding the results ofthe direct numerical simulation can be drawn: (i) first,the decay of the average polarization does not follow an7
FIG. 8. The distribution of dipole moments in direct numerical simulation in comparison to model predictions. Mea-surements are made after 1 [ µs ] of a step electric pulse (when stationarity is almost achieved). Distributions fol-low (left) a symmetric t -distribution in the transverse direction (along the y -axis) with model parameters ( ν, c, a, λ ) ≈ (5 . , . , .
108 [
A/mm ] , .
000 [
A/mm ]), while (right) in the direction parallel to the applied pulse (along the z -axis), weobserve a skewed t -distribution with model parameters ( ν, c, a, λ ) ≈ (7 . , . , .
164 [
A/mm ] , − .
864 [
A/mm ]). Note thatthe x -axes is the absolute value of the polarization. exponential decay; in fact it decays slightly slower thanan exponential function, (ii) second, under a constantapplied pulse, the variance increases initially but thendecreases to reach a plateau, and after switching off thepulse, variance exhibits an uptick before decaying to zero.The observed uptick in the variance is a unique featurethat is captured in our proposed model.We also solve for the model parameters in terms of theobserved distribution parameters, γ (cid:48) = (cid:115) ¯ γν − χ , (cid:15) = 1 (cid:113) a λ ,α (cid:48) = ¯ αγ (cid:48) (cid:15) ¯ γ − cγ (cid:48) √ − (cid:15) , u = − λγ (cid:48) (cid:15)α (cid:48) . The fitted values in figure 8 yield ( α (cid:48) , γ (cid:48) , (cid:15), u) =(2242 .
08 [ (cid:112) S / F] , . (cid:112) S / F] , . , , .
568 [A / mm ]),while ¯ γ = 1 . × [ S/F ] and ¯ α = 2 . × [ S/F ]. Fig-ure 9 depicts the comparison between our Fokker-Planckmodel and direct numerical simulations. We find that ourmodel perfectly captures the qualitative trends observedin the distribution of the polarizations, while it is in goodquantitative agreement. Also, it is important to developnumerical methods for the FPE in three spatial dimen-sions in order to faithfully compare these results, howeverthe current one dimensional analytic treatment providesvery encouraging results for the polarization componentparallel to the applied pulse. We will investigate numer-ical solutions to the full FPE in future works.
4. Fractional order evolution
In the case of Pearson diffusion, the statistical mo-ments up to order n < a − + 1 exist. In particular, for n ≥
2, the autocorrelation decays exponentially [109]: corr ( X s , X s + t ) = e − θt . In the current case, statistical moments exist up to order n < γγ (cid:48) . Therefore, this analysis suggests that the autocorrelationof dipole moments is given by an exponential decay when γ (cid:48) < γ. This result provides a critical threshold for the diver-sity measure γ (cid:48) c = √ γ that regulates the autocorrelationfunction, i.e. for γ (cid:48) < γ (cid:48) c , we predict anomalous relax-ation. Therefore, in our model, the observed anomalousrelaxation in cell aggregate electroporation is associatedwith the diversity in the cellular structural parametersthat is modeled through the parameter γ . In this case,the time derivative should be replaced with a fractionalorder derivative through ddt = τ α − d α dt α , i.e. τ is anarbitrary factor that has dimension of time. Here we8 FIG. 9. The evolution of average and variance of dipole moments in direct numerical simulations. We chose R = 7 [ µm ] and φ = 0 . × . i.e. a box of 4 [ mm ] on each side. The dottedmagenta line in the left figure illustrates the analytical solution (45). Letter N denotes numerical solution of FPE model whileA indicates analytical solution of the FPE model. assume τ = 1 and the governing equations read: d α µdt α = − ¯ γµ ( t ) − ¯ αu + γ (cid:48) µ ( t ) + (cid:15)γ (cid:48) α (cid:48) u, (43) d α σ dt α = − γ − γ (cid:48) ) σ ( t ) + γ (cid:48) µ ( t ) + 2 (cid:15)γ (cid:48) α (cid:48) uµ ( t ) + α (cid:48) u . (44)In order to preserve the type of initial conditions ap-propriate in classical phenomena, i.e. so that no extrainitial conditions be needed, we adopt Caputo fractionalderivative with m − < α ≤ m ( m is an integer number)[110, 111], which is defined by Ca D αt f ( t ) = 1Γ( m − α ) (cid:90) ta f ( m ) ( s )( t − s ) α − m +1 ds. Furthermore, defining fractional derivatives in the Ca-puto sense permits the application of the Laplace trans-form as a simple method of solution, see [112] for moredetails. In this case, the Laplace transform of Caputoderivatives reads: L [ C D αt f ( t )] = s α f ( s ) − m − (cid:88) k =0 f ( k ) (0 + ) s α − k − . Applying the Laplace transform to the set of equations(44) and assuming µ (0 + ) = 0 and µ (cid:48) (0 + ) = 0 yieldsthe transfer function H ( s ) (or impulse response) for theaverage polarization: µ ( s ) = H ( s ) u ( s ) , with H ( s ) = (cid:15)α (cid:48) γ (cid:48) − ¯ αs α + ¯ γ − γ (cid:48) . Then the impulse response function is given by: H ( t ) = ( (cid:15)α (cid:48) γ (cid:48) − ¯ α ) t α − E α,α [ − (¯ γ − γ (cid:48) ) t α ] , where E α,α is the Mittag-Leffler function that is generallydefined as: E α,β [ x ] = ∞ (cid:88) k =0 x k Γ( kα + β ) , α, β > . Note that for α = β = 1, it is equivalent to the expo-nential function. Therefore the general solution to theaverage polarization is given by: µ ( t ) (cid:15)α (cid:48) γ (cid:48) − ¯ α = (cid:90) t E α,α [ − (¯ γ − γ (cid:48) )( t − τ ) α ]( t − τ ) − α u ( τ ) dτ. Particularly, the step response waveform (in response to u ( t ) = u (0 + ) (cid:80) ∞ k =0 ( − k H ( t − t k ), where H ( t ) is the9 FIG. 10. Other time-domain properties based on our modelfor the case considered in figure 9. Immediately after switch-ing off the applied pulse a reverse current is observed.
Heaviside function) reads: µ ( t ) = u (0 + ) (cid:15)α (cid:48) γ (cid:48) − ¯ α ¯ γ − γ (cid:48) × ∞ (cid:88) k =0 ( − k (cid:18) − E α, [ − (¯ γ − γ (cid:48) )( t − t k ) α ] (cid:19) . (45)We found that a fractional order of α = 1 .
01 successfullydescribes the results obtained via direct numerical simu-lations, see figure 9. Alternatively, we could numericallyevaluate equation (44) using a finite difference numericalscheme [112, 113] which is basically to discretize Caputoderivative of order 0 < α < C D αt f ( t n +1 ) ≈ (∆ t ) − α Γ(2 − α ) n (cid:88) j =0 a j (cid:0) f n +1 − j − f n − j (cid:1) , where a j = ( j + 1) − α − j − α and f j = f ( t j ). Also, for1 < α < c.f. see equation 1.5of Sun and Wu (2006) [114]): C D αt f ( t n +1 ) ≈ (∆ t ) − α Γ(3 − α ) (cid:20) f n +1 − f n − b n − f (cid:48) (0)∆ t − n − (cid:88) j =1 ( b n − j − − b n − j )( f j − f j − ) (cid:21) , where b j = ( j +1) − α − j − α (see figure 9 for comparison). TABLE I. List of theoretical experiments based on the pro-posed framework. We consider a cell aggregate confined in acubic box with side length 1 [ mm ]. σ c [ S/m ] σ e [ S/m ] S L [ S/m ] φ α I 1 . . . . . . . . . . . × . . . . × . IV. PREDICTIONS & DISCUSSIONS
In this section we perform eight experiments based onthe proposed model with the specifications given in tablesI–II. In each experiment, we apply a Gaussian electricpulse given by: E ext ( t ) = E exp (cid:18) − t − t f / t f (cid:19) , where E = 40 [ kV /m ] for a duration of t f = 20 [ µs ]to resolve electric response at smaller frequencies; notethat the biggest frequency is determined by the max-imum time-step size of integration, which we limit to1 [ ns ], while the smallest frequency is inversely propor-tional to duration of integration. Then impedance canbe computed by equation (29).We performed numerical integrations of the integer or-der system of ordinary differential equations 41 with thepublicly available package Scipy [115] with adaptive timestepping. In particular, we used
LSODA algorithm, whichautomatically detects stiffness and switches between thenon-stiff Adam and stiff BDF integration methods [116].The fractional order ODEs 44 are solved by implementingthe discretization schemes discussed in subsection III C 4.The source code to solve the set of equations proposed inthis manuscript and to reproduce the results of this sec-tion can be found at https://github.com/pourion/CAEP.
A. Time-domain response
Figure 11 illustrates the evolution of the instantaneouseffective conductivity, current, resistance, and averagevalue of the membrane and the cytoplasm polarizationsduring the application of the external Gaussian pulse. Tomake these predictions, we considered an integer ordertime derivative, i.e. α = 1.At low membrane conductance (configurations I andII), we find that even though the magnitude of the cyto-plasm polarization is negligible with respect to the mem-brane polarization, over time the cytoplasm dipole mo-ment changes direction from anti-parallel to parallel withrespect to the external field, which leads to a slowly in-creasing resistance felt at the electrodes. However, in-creasing the membrane conductance (configurations III0 (a) Configuration I. (b) Configuration II.(c) Configuration III. (d) Configuration IV. FIG. 11. Experiments in the time domain. and IV) has the effect of increasing the cytoplasm polar-ization at the expense of reducing the membrane polar-ization while both dipole moments remain anti-parallelto the external field. Increasing the membrane conduc-tance enhances the overall current density and reducesthe overall electric resistance of the aggregate.We have shown how to compute the impedance di-rectly from the time-domain FPE, which paves the wayfor more detailed studies of cell aggregates with nonlinearmembrane processes such as the case of electroporation.
B. Impedance spectroscopy
The purpose of this section is to understand theimpedance as a function of frequency within cell ag-gregates. This analysis is important because it enablesthe resolution of the polarization processes and to relatethem to their relaxation timescales, cf. see the review byAsami (2002) [18]. For example impedance spectroscopyis widely used as a technique to characterize ionic con-1
TABLE II. List of theoretical experiments based on the pro-posed framework. We consider a cell aggregate confined in acubic box with side length 1 [ mm ]. In each case we chose 10different values for the varying parameter. σ c [ S/m ] σ e [ S/m ] S L [ S/m ] φ α V 0 . . . . , .
8] 1VI 1 . . , .
5] 1 . . . . . , . × ] 0 . . . . . . , . ductors, electroceramics, solid electrolytes, dielectric ma-terials such as polymers and glasses as well as fuel cellsand batteries [117–119].Figures 12–13 show the Bode and Cole diagrams cal-culated by equation (29) for 4 different configurationsand illustrate the effects caused by varying the volumefraction, the matrix conductivity, the membrane conduc-tance and the order of the fractional derivative.In figures 12(a)–12(b), we gradually increase the vol-ume fraction from φ = 1% to 80% that increases theimpedance at lower frequencies and reduces it at higherfrequencies. More importantly, increasing the volumefraction appears to amplify a low-frequency semi-circlein the Cole diagram originating from cell membranes.In figures 12(c)–12(d), we change the matrix conduc-tivity while keeping the other parameters fixed. We findthat when the matrix is more conductive than the cyto-plasm, the dielectric response of the cytoplasm lags be-hind that of the applied pulse. However for a cytoplasmmore conductive than the matrix, we find that the cy-toplam dielectric response leads the applied pulse. Thelatter behavior resembles the dielectric response of aninductive element that appears at high frequencies.In figures 13(a)–13(b), we gradually increase the mem-brane conductance and find that the semi-circle arc at lowfrequency gradually shrinks. The characteristic behaviorof the present model is that the membrane determinesthe low frequency arc while the cytoplasm determinesthe high frequency arc.In figures 13(c)–13(d), we vary the order of the frac-tional derivative. For 0 < α <
1, we observe a low fre-quency hook effect, where an apparently inductive loopappears at low frequencies, i.e. where the imaginary partof impedance becomes positive. In particular, we observethat our model predicts that, by increasing α towards 1,the low frequency hook gradually shrinks, and the lowfrequency semi-circle becomes depressed. Interestingly,Cole and Baker (1941) [120] reported an inductive re-sponse in their experiments with squid axons. Cole andBaker argued that inductive effects originated from themembrane of axons, which they modeled by an equiva-lent circuit composed of a resistor in series with an in-ductor that are connected in parallel with a capacitor.For a detailed discussion on the possible origins of in-ductive hooks we refer to Klotz (2019) [121]. In fact,low frequency inductive impedance is ubiquitously found in impedance spectroscopy experiments with various sys-tems such as Lithium ion batteries [122], proton exchangemembrane fuel cells [123], organic light emitting diodes(LEDs) [124], Perovskite solar cells [125], thin films onconductive substrates [126], and corrosion of Chromium[127]. It is well known that tissue impedance follows adepressed Cole (1940) equation, Z ( ω ) = R + R − R ∞ jω/ω ) α , where ω is the angular turnover frequency and α is a di-mensionless number between zero and one [26, 128, 129].It is generally established that it is the diversity of re-laxation timescales that is responsible for the observedanomalous electric response of tissue environments [130],which is the source of fractional order evolution in ourmodel as well. V. CONCLUSION
We have developed a theoretical framework based on adipole decomposition of cell polarization into two parts:the membrane polarization, and the cytoplasm polariza-tion. Based on this decomposition, we were able to eval-uate effective properties of the aggregate environmentsuch as effective conductivity and impedance. We alsoderived a time-domain governing Fokker-Planck equationthat explains distributions of cellular polarizations in dif-ferent volume fractions and at different frequencies. Weshowed that the effects of cell interactions can be easilyincluded in the model. Our theory is generally applica-ble to triphasic structures that are ubiquitously found innature, for example in modeling suspensions of biologicalcells and subcellular organella such as yeasts [131, 132],E. coli [133], synaptosomes [134], and mitochondria [135].The current work can be extended in several ways: • In plants and micro-organisms, cells are covered bya cell wall that adds another layer to the dielectricstructure; for details see Carstensen (1960) [136].Hanai et al. [137, 138] showed the number of inter-faces corresponds to the number of relaxations inthe dielectric response of a heterogeneous system,which could explain how diversity in the dielectricproperties of cells leads to anomalous relaxation.Therefore an extension of the current theory formultishell structures would be to develop N -phaseinterfacial polarization theories. • Coupling the bulk relaxation processes in tissue en-vironments, such as counterion polarization effects,with the interfacial polarization. In particular itwas argued [36] that counterion polarization effectscontribute to the observed anomalous relaxation;thus it will be useful to examine such influences onthe distribution of induced transmembrane poten-tials.2 (a) Configuration V - varying the volume fraction. (b) Configuration V - varying the volume fraction.(c) Configuration VI - varying the extra-cellular matrixconductivity. (d) Configuration VI - varying the extra-cellular matrixconductivity.
FIG. 12. Experiments in the frequency domain, in all figures warmer colors indicate higher values. Figure (a,b) show the effectof increasing volume fraction from φ = 0 .
01 to φ = 0 .
8. Figures (c,d) illustrate effects of increasing matrix conductivity in therange σ e = 0 . S/m ] to σ e = 1 . S/m ] while cytoplasm conductivity is fixed at σ c = 1 [ S/m ]. Red colors correspond to thecase of σ e > σ c while blue colors correspond to σ e < σ c . • Under strong electric fields, nonlinear cellular phe-nomena occur. A well known example is the mem-brane breakdown that occurs under transmem-brane potentials of about V ep = 0 . • Another interesting extension would be to considerthe effects of gap junctions on the induced trans-3 (a) Configuration VII - varying the membrane conductance. (b) Configuration VII - varying the membrane conductance.(c) Configuration VIII - varying the fractional order. (d) Configuration VIII - varying the fractional order.
FIG. 13. Experiments in the time domain, in all figures warmer colors correspond to higher values. Figures (a,b) capture theeffect of increasing the membrane conductance from S L = 1 . S/m ] (bluer colors) to S L = 1 . × [ S/m ] (redder colors).Figures (c,d) illustrate the effects of increasing the fractional order from α = 0 . α = 1 .
1. Note that the curves for α < α > membrane potentials. Gap junctions are electricalconnections between neighboring cells that providedirect pathways for ion transport in multicellularsystems. Gap junctions are key regulators for em-bryonic development due to their ability to regulatetransmembrane voltages; therefore understandingtheir interplay with an external electric stimula- tion poses new opportunities to control embryonicdevelopment and a new pathway to understand andcontrol patterning in biological organisms.4
ACKNOWLEDGMENTS
Part of this research was funded by ARO W911NF-16-1-0136. The direct numerical simulations performedin this work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported byNational Science Foundation grant number ACI-1053575and ASC-150002. The authors acknowledge the TexasAdvanced Computing Center (TACC) at The Universityof Texas at Austin for providing HPC and visualizationresources that have contributed to the research resultsreported within this paper. [1] U. Zimmermann, Electric field-mediated fusion andrelated electrical phenomena, Biochimica et Biophys-ica Acta (BBA)-Reviews on Biomembranes , 227(1982).[2] C. A. Jordan, E. Neumann, and A. E. Sowers,
Electropo-ration and electrofusion in cell biology (Springer Science& Business Media, 2013).[3] W. Arnold and U. Zimmermann, Electric field-inducedfusion and rotation of cells, Biological membranes , 389(1984).[4] G. Fuhr, R. Glaser, and R. Hagedorn, Rotation ofdielectrics in a rotating electric high-frequency field.model experiments and theoretical explanation of therotation effect of living cells, Biophysical journal ,395 (1986).[5] H. A. Pohl and J. S. Crane, Dielectrophoresis of cells,Biophysical journal , 711 (1971).[6] F. A. Sauer, Forces on suspended particles in the elec-tromagnetic field, in Coherent excitations in biologicalsystems (Springer, 1983) pp. 134–144.[7] F. F. Becker, X.-B. Wang, Y. Huang, R. Pethig, J. Vyk-oukal, and P. Gascoyne, Separation of human breastcancer cells from blood by differential dielectric affin-ity, Proceedings of the National Academy of Sciences , 860 (1995).[8] P. R. Gascoyne and J. Vykoukal, Particle separation bydielectrophoresis, Electrophoresis , 1973 (2002).[9] E. Neumann and K. Rosenheck, Permeability changesinduced by electric impulses in vesicular membranes,The Journal of membrane biology , 279 (1972).[10] K. Kaler and T. Jones, Dielectrophoretic spectra of sin-gle cells determined by feedback-controlled levitation.,Biophysical journal , 173 (1990).[11] G. Bryant and J. Wolfe, Electromechanical stresses pro-duced in the plasma membranes of suspended cells byapplied electric fields, The Journal of membrane biology , 129 (1987).[12] G. H. Markx and C. L. Davey, The dielectric propertiesof biological cells at radiofrequencies: applications inbiotechnology, Enzyme and Microbial Technology ,161 (1999).[13] M. Levin, G. Pezzulo, and J. M. Finkelstein, Endoge-nous bioelectric signaling networks: exploiting voltagegradients for control of growth and form, Annual reviewof biomedical engineering , 353 (2017).[14] W. Roux, Uber die* morphologische˜ polarisation yoneiern und embryonen durch den elektrischen strom,scwie fiber die wirkung des elektrischen stromes auf dierichtung der ersten teilung des eies, Sitzungsber. der kkAkad. d. Wiss. Wien math. nat. K , 27.[15] K. B. Hotary and K. R. Robinson, Evidence of a rolefor endogenous electrical fields in chick embryo develop- ment, Development , 985 (1992).[16] C. E. Pullar, The physiology of bioelectricity in devel-opment, tissue regeneration and cancer (CRC Press,2016).[17] B. Reid and M. Zhao, The electrical response to injury:molecular mechanisms and wound healing, Advances inwound care , 184 (2014).[18] K. Asami, Characterization of heterogeneous systemsby dielectric spectroscopy, Progress in polymer science , 1617 (2002).[19] P. Mistani, S. Pakravan, and F. Gibou, A fractionalstochastic theory for nonlinear electroporation of cellaggregates, in preparation.[20] J. Cervera, M. Levin, and S. Mafe, Bioelectrical cou-pling of single-cell states in multicellular systems, TheJournal of Physical Chemistry Letters , 3234 (2020).[21] H. P. Schwan, Electrical properties of tissue and cell sus-pensions, in Advances in biological and medical physics ,Vol. 5 (Elsevier, 1957) pp. 147–209.[22] M. Stuchly and S. Stuchly, Dielectric properties ofbiological substances?tabulated, Journal of Microwavepower , 19 (1980).[23] R. Pethig, Dielectric properties of biological materials:Biophysical and medical applications, IEEE Transac-tions on Electrical Insulation , 453 (1984).[24] R. Pethig and D. B. Kell, The passive electrical proper-ties of biological systems: their significance in physiol-ogy, biophysics and biotechnology, Physics in Medicine& Biology , 933 (1987).[25] K. R. Foster, H. P. Schwan, et al. , Dielectric propertiesof tissues, CRC handbook of biological effects of elec-tromagnetic fields , 27.[26] E. McAdams and J. Jossinet, Tissue impedance: ahistorical overview, Physiological measurement , A1(1995).[27] C. Gabriel, S. Gabriel, and y. E. Corthout, The dielec-tric properties of biological tissues: I. literature survey,Physics in medicine & biology , 2231 (1996).[28] W. Kuang and S. Nelson, Low-frequency dielectric prop-erties of biological tissues: a review with some new in-sights, (1998).[29] T. Kotnik, L. Rems, M. Tarek, and D. Miklavˇciˇc,Membrane electroporation and electropermeabilization:mechanisms and models, Annual review of biophysics , 63 (2019).[30] H. Schwan, G. Schwarz, J. Maczuk, and H. Pauly, Onthe low-frequency dielectric dispersion of colloidal par-ticles in electrolyte solution1, The Journal of PhysicalChemistry , 2626 (1962).[31] G. Schwarz, A theory of the low-frequency dielectric dis-persion of colloidal particles in electrolyte solution1, 2,The Journal of Physical Chemistry , 2636 (1962). [32] C. W. Einolf Jr and E. L. Carstensen, Low-frequency di-electric dispersion in suspensions of ion-exchange resins,The Journal of Physical Chemistry , 1091 (1971).[33] C. W. Einolf Jr and E. L. Carstensen, Passive electri-cal properties of microorganisms: V. low-frequency di-electric dispersion of bacteria, Biophysical journal , 8(1973).[34] S. S. Dukhin, V. N. Shilov, and J. Bikerman, Dielectricphenomena and double layer in disperse systems andpolyelectrolytes, Journal of the Electrochemical Society , 154C (1974).[35] M. Mandel and T. Odijk, Dielectric properties of poly-electrolyte solutions, Annual Review of Physical Chem-istry , 75 (1984).[36] C. Grosse and K. R. Foster, Permittivity of a suspen-sion of charged spherical particles in electrolyte solution,Journal of Physical Chemistry , 3073 (1987).[37] J. C. Maxwell, A treatise on electricity and magnetism ,Vol. 1 (Oxford: Clarendon Press, 1873).[38] K. W. Wagner, Erkl¨arung der dielektrischen nach-wirkungsvorg¨ange auf grund maxwellscher vorstellun-gen, Archiv f¨ur Elektrotechnik , 371 (1914).[39] H. Fricke, A mathematical treatment of the electric con-ductivity and capacity of disperse systems ii. the capac-ity of a suspension of conducting spheroids surroundedby a non-conducting membrane for a current of low fre-quency, Physical Review , 678 (1925).[40] H. Fricke, The electric permittivity of a dilute suspen-sion of membrane-covered ellipsoids, Journal of AppliedPhysics , 644 (1953).[41] T. Hanai, Theory of the dielectric dispersion due to theinterfacial polarization and its application to emulsions,Kolloid-Zeitschrift , 23 (1960).[42] T. Hanai, K. Asami, and N. Koizumi, Dielectric theoryof concentrated suspensions of shell-spheres in partic-ular reference to the analysis of biological cell suspen-sions, (1979).[43] H. Zhang, K. Sekine, T. Hanai, and N. Koizumi, Dielec-tric observations on polystyrene microcapsules and thetheoretical analysis with reference to interfacial polar-ization, Colloid and Polymer Science , 381 (1983).[44] H. Zhang, K. Sekine, T. Hanai, and N. Koizumi, Dielec-tric approach to polystyrene microcapsule analysis andthe application to the capsule permeability to potassiumchloride, Colloid and polymer science , 513 (1984).[45] T. Hanai, K. Zhao, K. Asaka, and K. Asami, Theo-retical approach and the practice to the evaluation ofstructural parameters characterizing concentration po-larization alongside ion-exchange membranes by meansof dielectric measurement, Colloid and Polymer Science , 766 (1993).[46] P. J. W. Debye, Polar molecules (Chemical CatalogCompany, Incorporated, 1929).[47] J. Miles Jr and H. Robertson, The dielectric behavior ofcolloidal particles with an electric double-layer, PhysicalReview , 583 (1932).[48] T. Kotnik, F. Bobanoviˇc, and D. Miklavˇciˇc, Sensitivityof transmembrane voltage induced by applied electricfields – a theoretical analysis, Bioelectrochemistry andbioenergetics , 285 (1997).[49] D. Gross, Electromobile surface charge alters membranepotential changes induced by applied electric fields, Bio-physical journal , 879 (1988).[50] E. Neumann, The relaxation hysteresis of membrane electroporation, in Electroporation and Electrofusion inCell Biology , edited by E. Neumann, A. E. Sowers, andC. A. Jordan (Springer US, Boston, MA, 1989) pp. 61–82.[51] S. Ho and G. Mittal, Electroporation of cell mem-branes: a review, Critical reviews in biotechnology ,349 (1996).[52] T. Geng and C. Lu, Microfluidic electroporation for cel-lular analysis and delivery, Lab on a Chip , 3803(2013).[53] T. Kotnik and D. Miklavˇciˇc, Analytical description oftransmembrane voltage induced by electric fields onspheroidal cells, Biophysical Journal , 670 (2000).[54] J. Teissie and M.-P. Rols, An experimental evaluationof the critical potential difference inducing cell mem-brane electropermeabilization, Biophysical journal ,409 (1993).[55] M. Legu`ebe, A. Silve, L. Mir, and C. Poignard,Conducting and permeable states of cell mem-brane submitted to high voltage pulses: Math-ematical and numerical studies validated bythe experiments, Journal of Theoretical Biologyhttps://doi.org/10.1016/j.jtbi.2014.06.027 (2014).[56] A. Guittet, M. Lepilliez, S. Tanguy, and F. Gibou,Solving elliptic problems with discontinuities on irregu-lar domains – the voronoi interface method, Journal ofComputational Physics , 747 (2015).[57] R. P. Fedkiw, T. Aslam, B. Merriman, S. Osher, et al. ,A non-oscillatory eulerian approach to interfaces in mul-timaterial flows (the ghost fluid method), (1999).[58] A. Guittet, C. Poignard, and F. Gibou, A voronoi inter-face approach to cell aggregate electropermeabilization,Journal of Computational Physics , 143 (2017).[59] P. Mistani, A. Guittet, C. Poignard, and F. Gibou, Aparallel voronoi-based approach for mesoscale simula-tions of cell aggregate electropermeabilization, Journalof Computational Physics , 48 (2019).[60] S. Balay, S. Abhyankar, M. F. Adams, J. Brown,P. Brune, K. Buschelman, L. Dalcin, A. Dener, V. Ei-jkhout, W. D. Gropp, D. Karpeyev, D. Kaushik,M. G. Knepley, D. A. May, L. C. McInnes, R. T.Mills, T. Munson, K. Rupp, P. Sanan, B. F. Smith,S. Zampini, H. Zhang, and H. Zhang, PETSc Web page, (2019).[61] S. Balay, S. Abhyankar, M. F. Adams, J. Brown,P. Brune, K. Buschelman, L. Dalcin, A. Dener, V. Ei-jkhout, W. D. Gropp, D. Karpeyev, D. Kaushik,M. G. Knepley, D. A. May, L. C. McInnes, R. T.Mills, T. Munson, K. Rupp, P. Sanan, B. F. Smith,S. Zampini, H. Zhang, and H. Zhang, PETSc UsersManual , Tech. Rep. ANL-95/11 - Revision 3.13 (Ar-gonne National Laboratory, 2020).[62] S. Balay, W. D. Gropp, L. C. McInnes, and B. F. Smith,Efficient management of parallelism in object orientednumerical software libraries, in
Modern Software Toolsin Scientific Computing , edited by E. Arge, A. M. Bru-aset, and H. P. Langtangen (Birkh¨auser Press, 1997) pp.163–202.[63] H. A. Van der Vorst, Bi-cgstab: A fast and smoothlyconverging variant of bi-cg for the solution of nonsym-metric linear systems, SIAM Journal on scientific andStatistical Computing , 631 (1992).[64] R. D. Falgout and U. M. Yang, hypre: A library ofhigh performance preconditioners, in International Con- ference on Computational Science (Springer, 2002) pp.632–641.[65] C. Burstedde, L. C. Wilcox, and O. Ghattas, p4est :Scalable algorithms for parallel adaptive mesh refine-ment on forests of octrees, SIAM Journal on ScientificComputing , 1103 (2011).[66] C. Rycroft, Voro++: A three-dimensional Voronoicell library in C++ , Tech. Rep. (Lawrence BerkeleyNational Lab.(LBNL), Berkeley, CA (United States),2009).[67] T. C. Choy,
Effective medium theory: principles andapplications , Vol. 165 (Oxford University Press, 2015).[68] A. D. Fokker, Die mittlere energie rotierender elek-trischer dipole im strahlungsfeld, Annalen der Physik , 810 (1914).[69] M. Planck, Uber einen satz der statistischen dynamikund seine erweiterung in der quantentheorie, , 324(1917).[70] D. J. Jeffrey, Conduction through a random suspensionof spheres, Proceedings of the Royal Society of London.A. Mathematical and Physical Sciences , 355 (1973).[71] G. Batchelor, Sedimentation in a dilute dispersion ofspheres, Journal of fluid mechanics , 245 (1972).[72] G. Batchelor and J.-T. Green, The hydrodynamic in-teraction of two small freely-moving spheres in a linearflow field, Journal of Fluid Mechanics , 375 (1972).[73] P. C. Franzone and G. Savar´e, Degenerate evolution sys-tems modeling the cardiac electric field at micro-andmacroscopic level, in Evolution equations, semigroupsand functional analysis (Springer, 2002) pp. 49–78.[74] A. Collin and S. Imperiale, Mathematical analysisand 2-scale convergence of a heterogeneous microscopicbidomain model, Mathematical Models and Methods inApplied Sciences , 979 (2018).[75] J. Hinch, A perspective of batchelor’s research in micro-hydrodynamics, Journal of Fluid Mechanics , 8(2010).[76] G. Batchelor, Transport properties of two-phase mate-rials with random structure, Annual Review of FluidMechanics , 227 (1974).[77] R. Bonnecaze and J. Brady, A method for determin-ing the effective conductivity of dispersions of particles,Proceedings of the Royal Society of London. Series A:Mathematical and Physical Sciences , 285 (1990).[78] H. Risken, The fokker-planck equation: methods of so-lution and applications, (1984).[79] W. Coffey and Y. P. Kalmykov, The Langevin equa-tion: with applications to stochastic problems in physics,chemistry and electrical engineering , Vol. 27 (World Sci-entific, 2012).[80] J. D. Jackson, Classical electrodynamics (1999).[81] D. B. Geselowitz, On bioelectric potentials in an inho-mogeneous volume conductor, Biophysical journal , 1(1967).[82] W. R. Smythe, Static and dynamic electricity , 3rded., International series in pure and applied physics(McGraw-Hill, New York, 1967).[83] Y. Chiew and E. Glandt, The effect of structure on theconductivity of a dispersion, Journal of colloid and in-terface science , 90 (1983).[84] D. Hasselman and L. F. Johnson, Effective thermal con-ductivity of composites with interfacial thermal bar-rier resistance, Journal of composite materials , 508(1987). [85] Y. C. Chiew and E. D. Glandt, Effective conductivityof dispersions: the effect of resistance at the particlesurfaces, Chemical engineering science , 2677 (1987).[86] R. O’brien, A method for the calculation of the effec-tive transport properties of suspensions of interactingparticles, Journal of Fluid Mechanics , 17 (1979).[87] J. Dahler and L. Scriven, Theory of structured continuai. general consideration of angular momentum and po-larization, Proceedings of the Royal Society of London.Series A. Mathematical and Physical Sciences , 504(1963).[88] G. Russakoff, A derivation of the macroscopic maxwellequations, American Journal of Physics , 1188 (1970).[89] J. C. Maxwell, A treatise on electricity and magnetism ,Vol. 1 (Clarendon press, 1881).[90] D. Ross, The potential due to two point charges each atthe centre of a spherical cavity and embedded in a di-electric medium, Australian Journal of Physics , 817(1968).[91] S.-Y. Lu and J.-L. Song, Effective conductivity of com-posites with spherical inclusions: Effect of coating anddetachment, Journal of applied physics , 609 (1996).[92] Z. Hashin and S. Shtrikman, A variational approach tothe theory of the effective magnetic permeability of mul-tiphase materials, Journal of applied Physics , 3125(1962).[93] P. Langevin, Sur la th´eorie du mouvement brownien,Compt. Rendus , 530 (1908).[94] E. Nelson, Dynamical theories of Brownian motion ,Vol. 3 (Princeton university press, 1967).[95] H. Takayasu, A.-H. Sato, and M. Takayasu, Stableinfinite variance fluctuations in randomly amplifiedlangevin systems, Physical Review Letters , 966(1997).[96] G. Steinbrecher and B. Weyssow, Extreme anoma-lous particle transport at the plasma edge, Univ. ofCraiova/Univ Libre de Bruxells working paper (2007).[97] A. Schenzle and H. Brand, Multiplicative stochastic pro-cesses in statistical physics, Physical Review A , 1628(1979).[98] R. M. May, S. A. Levin, and G. Sugihara, Ecology forbankers, Nature , 893 (2008).[99] W. T. Shaw and M. Schofield, A model of returns forthe post-credit-crunch reality: hybrid brownian mo-tion with price feedback, Quantitative Finance , 975(2015).[100] J. L. Forman and M. Sørensen, The pearson diffusions:A class of statistically tractable diffusion processes,Scandinavian Journal of Statistics , 438 (2008).[101] H. Risken, Fokker-planck equation, in The Fokker-Planck Equation: Methods of Solution and Applications (Springer Berlin Heidelberg, Berlin, Heidelberg, 1996)pp. 63–95.[102] H. Hasegawa, Stationary and dynamical properties ofinformation entropies in nonextensive systems, PhysicalReview E , 031133 (2008).[103] H. Hasegawa, A moment approach to analytic time-dependent solutions of the fokker-planck equationwith additive and multiplicative noise, arXiv preprintarXiv:0810.0839 (2008).[104] R. E. Mortensen, Mathematical problems of modelingstochastic nonlinear dynamic systems, Journal of Sta-tistical Physics , 271 (1969).[105] K. Pearson, Tables for statisticians and biometricians (University press, 1914).[106] S. M. Iacus, Simulation and inference for stochastic dif-ferential equations: with R examples (Springer Science& Business Media, 2009).[107] K. Pearson, X. contributions to the mathematical the-ory of evolution.?ii. skew variation in homogeneous ma-terial, Philosophical Transactions of the Royal Societyof London.(A.) , 343 (1895).[108] J. Heinrich, A guide to the pearson type iv distribution,(2004).[109] B. M. Bibby, I. M. Skovgaard, M. Sørensen, et al. ,Diffusion-type models with given marginal distributionand autocorrelation function, Bernoulli , 191 (2005).[110] M. Caputo, Linear models of dissipation whose q is al-most frequency independent?ii, Geophysical Journal In-ternational , 529 (1967).[111] R. Gorenflo and F. Mainardi, Fractional calculus: inte-gral and differential equations of fractional order, arXivpreprint arXiv:0805.3823 (2008).[112] M. A. Matlob and Y. Jamali, The concepts and applica-tions of fractional order differential calculus in modelingof viscoelastic systems: A primer, Critical Reviews? inBiomedical Engineering (2019).[113] C. Li and F. Zeng, Finite difference methods for frac-tional differential equations, International Journal of Bi-furcation and Chaos , 1230014 (2012).[114] Z.-z. Sun and X. Wu, A fully discrete difference schemefor a diffusion-wave system, Applied Numerical Mathe-matics , 193 (2006).[115] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haber-land, T. Reddy, D. Cournapeau, E. Burovski, P. Pe-terson, W. Weckesser, J. Bright, S. J. van der Walt,M. Brett, J. Wilson, K. Jarrod Millman, N. Mayorov,A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. Carey,˙I. Polat, Y. Feng, E. W. Moore, J. Vand erPlas, D. Lax-alde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quin-tero, C. R. Harris, A. M. Archibald, A. H. Ribeiro,F. Pedregosa, P. van Mulbregt, and S. . . Contributors,SciPy 1.0: Fundamental Algorithms for Scientific Com-puting in Python, Nature Methods , 261 (2020).[116] L. Petzold, Automatic selection of methods for solvingstiff and nonstiff systems of ordinary differential equa-tions, SIAM journal on scientific and statistical comput-ing , 136 (1983).[117] J. R. Macdonald and E. Barsoukov, Impedance spec-troscopy: theory, experiment, and applications, History , 1 (2005).[118] M. E. Orazem and B. Tribollet, Electrochemicalimpedance spectroscopy, New Jersey (2008).[119] N. Bonanos, P. Pissis, and J. R. Macdonald, Impedancespectroscopy of dielectrics and electronic conductors,Characterization of materials , 1 (2002).[120] K. S. Cole and R. F. Baker, Longitudinal impedance ofthe squid giant axon, The Journal of general physiology , 771 (1941).[121] D. Klotz, Negative capacitance or inductive loop?–a general assessment of a common low frequencyimpedance feature, Electrochemistry Communications , 58 (2019). [122] H. Brandst¨atter, I. Hanzu, and M. Wilkening, Myth andreality about the origin of inductive loops in impedancespectra of lithium-ion electrodes?a critical experimentalapproach, Electrochimica acta , 218 (2016).[123] S. K. Roy, M. E. Orazem, and B. Tribollet, Interpreta-tion of low-frequency inductive loops in pem fuel cells,Journal of The Electrochemical Society , B1378(2007).[124] J. Bisquert, G. Garcia-Belmonte, ´A. Pitarch, and H. J.Bolink, Negative capacitance caused by electron injec-tion through interfacial states in organic light-emittingdiodes, Chemical Physics Letters , 184 (2006).[125] E. Ghahremanirad, A. Bou, S. Olyaee, and J. Bisquert,Inductive loop in the impedance response of perovskitesolar cells explained by surface polarization model, Thejournal of physical chemistry letters , 1402 (2017).[126] S. Taibl, G. Fafilek, and J. Fleig, Impedance spectraof fe-doped srtio 3 thin films upon bias voltage: induc-tive loops as a trace of ion motion, Nanoscale , 13954(2016).[127] J. Dobbelaar and J. De Wit, Impedance measurementsand analysis of the corrosion of chromium, Journal ofThe Electrochemical Society , 2038 (1990).[128] K. S. Cole and R. H. Cole, Dispersion and absorptionin dielectrics i. alternating current characteristics, TheJournal of chemical physics , 341 (1941).[129] S. Smye, C. Evans, M. Robinson, and B. Sleeman, Mod-elling the electrical properties of tissue as a porousmedium, Physics in Medicine & Biology , 7007(2007).[130] R. Zorn, Logarithmic moments of relaxation time dis-tributions, The Journal of chemical physics , 3204(2002).[131] Y. SUGIURA, S. KOGA, and H. AKABORI, Dielectricbehavior of yeast cells in suspension, The Journal ofGeneral and Applied Microbiology , 163 (1964).[132] K. Asami, T. Hanai, and N. Koizumi, Dielectric prop-erties of yeast cells, The Journal of membrane biology , 169 (1976).[133] H. Fricke, H. P. Schwan, K. Li, and V. Bryson, A dielec-tric study of the low-conductance surface membrane ine. coli, Nature , 134 (1956).[134] A. Irimajiri, T. Hanai, and A. Inouye, Dielectric prop-erties of synaptosomes isolated from rat brain cortex,Biophysics of structure and mechanism , 273 (1975).[135] H. Pauly, L. Packer, and H. Schwan, Electrical proper-ties of mitochondrial membranes, The Journal of bio-physical and biochemical cytology , 589 (1960).[136] E. Carstensen and R. Marquis, Passive electrical prop-erties of microorganisms: Iii. conductivity of isolatedbacterial cell walls, Biophysical journal , 536 (1968).[137] T. Hanai and K. Sekine, Theory of dielectric relaxationsdue to the interfacial polarization for two-componentsuspensions of spheres, Colloid and Polymer Science , 888 (1986).[138] T. Hanai, H. Zhang, K. Sekine, K. Asaka, and K. Asami,The number of interfaces and the associated dielectricrelaxations in heterogeneous systems, Ferroelectrics86