You are viewing a javascript disabled version of the site. Please enable Javascript for this site to function properly.

# Spectroscopic conductivity imaging of a cell culture

#### Abstract

In this paper, we present a simplified electrical model for tissue culture. We derive a mathematical structure for overall electrical properties of the culture and study their dependence on the frequency of the current. We introduce a method for recovering the microscopic properties of the cell culture from the spectral measurements of the effective conductivity. Numerical examples are provided to illustrate the performance of our approach.

## 1.Introduction

Cell culture production processes, such as those from stem cell therapy, must be monitored and controlled to meet strict functional requirements. For example, a cell culture of cartilage, designed to replace that in the knee, must be organized in a specific way.

Hyaline cartilage is located on the joint surface and play an important role in body movement. In normal articular cartilage, there is a depth-dependent stratified structure known as zonal organization. As a simplified model, cartilage comprises three different layers [10]: a superficial zone in outer 10%, a middle zone that is 50% of the height, and a deep zone consisting in the inner 40%. At the microscopic level, cartilage tissue is composed of cells, collagen fibers, and glycosaminoglycans (GAGs). The concentration and organization of each microstructure differs among the three layers. In the superficial zone, cells are anisotropic and horizontally aligned, collagen orientation is also horizontal and GAGs have a lower concentration than in the other layers. In the middle zone, there are fewer cells and they are isotropic, collagen is randomly oriented and there is a medium concentration of GAGs. In the deep zone, cells are isotropic, cell density is higher than in the middle zone, collagen is vertically aligned and there is a high GAG density. As these parameters all contribute to the function of collagen in the knee, and must be replicated in the cell culture.

It is important that the method for monitoring cell cultures is non-destructive. Destructive methods require hundreds of samples to be cultured for a single functional tissue, and for the samples to be monitored multiple times during maturation. Here, we propose a microscopic electrical impedance tomography (micro-EIT) method for monitoring cell cultures that exploits the distinctive dielectric properties of cells and other microstructures. In this method, electrodes inject a current into the medium at different frequencies and the corresponding dielectric potentials are recorded, thus enabling reconstruction of the microscopic parameters of the medium. The parameters of interest are cell density, collagen orientation, and GAG density, as well as the orientation and shape of cells.

EIT uses a low-frequency current (below 500 kHz) to visualize the internal impedance distribution of a conducting domain such as a tissue sample or the human body. Recent studies measured electrical conductivity values and anisotropy ratios of engineered cartilage to distinguish extracellular matrix samples containing differing amounts of collagen and GAGs. During chondrogenesis over a six-week period, these measurements could distinguish the stages of the process and provide information regarding the internal depth-dependent structure.

In this work, we provide a mathematical framework for determining the microscopic properties of the cell culture from spectral measurements of the effective conductivity. For simplicity, we consider a microstructure comprising two components in a background medium. One of the components has a frequency dependent on the material parameters arising from the cell membrane structure, while the other has constant conductivity and permittivity over the frequency range. First, we derive in Theorem 2 the overall electrical properties of the culture, which depend on the volume fraction of each component and associated membrane polarization tensors defined by (10) and (11). Then, we show that the spectral measurements of the overall electrical properties of the culture can be used to determine the volume fraction of each component and the anisotropy ratio of the first component. For doing so, we study the dependence of the membrane polarization tensors on the operating frequency and use the spectral theorem to recover in Proposition 9 from the measurement of the effective conductivity on a range of frequencies the coefficients of its expansion with respect to the frequency. Proposition 9 also provides the anisotropy ratio of the cell culture.

This paper is organized as follows. In Section 2, we present a simplified model of the tissue culture. In Section 3, we derive an equivalent effective conductivity for the solution at the macroscopic scale. In Section 4, we present a method based on spectral measurements, in which microscopic properties are measured from the effective conductivity. This process is known as inverse homogenization or dehomogenization [5,14]. Finally, we provide some numerical examples to illustrate our main findings.

## 2.The direct problem

In this section, we propose a simple electrical model for the tissue and derive an effective conductivity using periodic homogenization.

### 2.1.Problem setting

We consider the domain of interest – the cell culture – to be described by a domain ΩR3. We assume that Ω=D×(0,1) where DR2 denotes a floor of the culture medium; see Fig. 1. Following [11,17,18], we describe the conductivity of the medium by a scalar field

σω,ϵ(x)=σω(x,xϵ),
where ω denotes the angular frequency of the injected current, and ϵ>0 is a small parameter representing the microscopic scale of the medium; σ is 1-periodic in every direction in the second variable. Let us consider the following unit domain:
Y=(12;12)d.
For a fixed x, σ(x,xϵ) describes the conductivity in a single cartilage tissue with cell size ϵ at a location xΩ. To have a complete model of the tissue, σ must describe the conductivity of both cells and of the other inclusions, i.e., collagen and GAGs. The biological fluid conductivity is noted k0 and is assumed to be frequency independent. The cells are made of biological fluid enclosed in a very thin and very resistive membrane [2] of thickness ϵδ for some small parameter δ>0. The conductivity of the membrane is frequency dependent and is noted km(ω). The cell shape varies slowly with the parameter xΩ compared to the microscale ϵ. The other inclusions are described by some frequency independent conductivity function ki(x,xϵ). Let
ψ:Ω×Rd:R
be a C1(Ω×Rd) function, 1-periodic in every direction with respect to the second variable. We assume that the function ψ is the level set function for the membrane boundary given by Ωϵ+={x:ψ(x,xϵ)>δ} (resp. Ωϵ={x:ψ(x,xϵ)<δ}). We also assume that the support of ki(x,y) is strictly included in {(x,y):ψ(x,y)>δ}. We can now describe the conductivity σω, which is schematically represented at a fixed x in Fig. 2:
(1)σω(x,y)=k0+ki(x,y)if ψ(x,y)>δ,k0if ψ(x,y)<δ,km(ω)else.

##### Fig. 1.

Organization of the cells in the cartilage tissue.

##### Fig. 2.

Typical values of σω on Y.

Now that we have an expression for the conductivity in the medium, as commonly accepted in EIT, we use the quasistatic approximation for the electrical potential. For an input current g(x)sin(ωt) on the boundary Y, with Ωg=0, the real part of the corresponding time-harmonic potential, denoted by uω,ϵ, satisfies the following problem approximately:

(2)·σω,ϵuω,ϵ=0in Ω,σω,ϵuω,ϵ·ν=gon Ω.
where ν is the outer normal vector on Ω. Here, we impose the normalization Ωϵuω,ϵ=0.

##### Remark 1.

Let us briefly explain how the expression of σω in (1) is derived. We should note that the frequency dependent behaviors of σω,ϵ in (2) are attributed to thin cell membranes. Imagine that we inject an oscillating current at the angular frequency ω into the cube Y. Then, the resulting time-harmonic potential w=u+iv in Y is governed by

·((σ(y)+iωσ(y))w(y))=0for yY,
where σ denotes the conductivity distribution and σ is the permittivity distribution in Y. In [9], it was shown that, under some conditions on the membrane, the real part u approximately satisfies
(3)·(|σ+iωσ|2σu)=0in Y.
Since σσ outside the membrane, we have
|σ+iωσ|2σσoutside the membrane.
Hence, it is reasonable to assume that the conductivity outside the membrane, as a coefficient of the elliptic PDE (3), does not change with frequency. On the other hand, since σ on the membrane is very small, the effect of σ is not negligible. Hence, the conductivity, km, on the membrane changes with frequency as follows:
|σ+iωσ|2σ=σ+ω2σσon the membrane.

### 2.2.Homogenization of the tissue

We are now interested in getting rid of the microscale oscillations of σω,ϵ, since boundary measurements will only allow us to image macroscale variations of the conductivity. To this end, we proceed to the homogenization of equation (2). Assume that k0+ki is bounded from below and from above:

0<σ_k0+kiσ.
From [2], we have two-scale convergence [1,11,13] of uω,ϵ to uω, which is a solution to
(4)·σωuω=0in Ω,σωuω·ν=gon Ω,Ωuω=0,
for an input current g(x)sin(ωt) on the boundary Ω. Here, σω is called the effective conductivity which can be represented by [2]
σω(x)ep·eq=Yσω(x,y)(yp+vp(y))·eqdy,p,q{1,,d}(5)=k0(δp,q+Yvpνyqds(y)),
where ep:=(0,,1,,0) with 0 components except pth component 1, and vp is the solution to the following equation on Y for p=1,,d:
(6)·(σω(x,y)(vp(y)+yp))=0for yY,vp1-periodic,Y(vp(y)+yp)dy=0.
As δ0, vp can be approximated [8,15,16] by the solution of the following equation, where β(ω)=δkm(ω):
(7)·(σω(x,y)(vp(y)+yp))=0for yYC,k0ν(vp+(y)+yp)=k0ν(vp(y)+yp)for yC,vp+(y)vp(y)=β(ω)k0ν(vp+(y)+yp)for yC,vp1-periodic,Yvp(y)+ydy=0.
Here, C denotes the membrane of the cell C and βk0 is the effective thickness of the membrane.

## 3.Imaging the microstructure from effective conductivity measurements

In this section, we do not care about the space dependence of σω, and will therefore drop it. We will thus assume that σω is constant equal to some matrix in Md(C):={mCd×d:mi,j=mj,i for i,j=1,2,,d}. We will show what kind of information on the microstructure we can recover from the knowledge of σω in a range of frequencies ω(ω1,ω2). First, in Section 3.1, we will obtain a simple representation of the effective conductivity in the dilute case, where the volume fraction of both cells and other inclusions is small compared to the volume of biological fluid. Then, in the following sections we will use this representation and will show how to recover information about the microstructure using the spectral measure.

### 3.1.Effective conductivity in the dilute case

Here, we consider some reference cell C0 and some reference inclusion B0 with there C2 boundaries C0 and B0. We assume that C=xC+ρCC0 and β(ω)=ρCβ0(ω) for some reference β0(ω) and let B=xB+ρBB0, where xC and xB respectively indicate the locations of the cell and inclusion and ρC and ρB their characteristic sizes. We assume that the conductivity ki of the inclusion is given by

ki(y)=(k0k1)χB(y),
where χB denotes the characteristic function of B.

The effective conductivity is therefore expressed as

σωep·eq=Yσ(y)(yp+vp(y))·eqdy,p,q{1,,d},
where, for p{1,,d},
(8)·(k0(vp(y)+yp))=0in Y(BC),·(k1(vp(y)+yp))=0in B,k0ν(vp+(y)+yp)=k0ν(vp(y)+yp)on C,vp+vp=β(ω)k0ν(vp+(y)+yp)on C,vp+vp=0on B,k0ν(vp+(y)+yp)=k1ν(vp(y)+yp)on B,vpperiodic,Y(vp(y)+y)dy=0.
From now on, I denotes the inclusion map H1/2(C)H1/2(C), where H1/2 and H1/2 are the Sobolev spaces of order 1/2 and 1/2 on C. We will now proceed to prove the following result.

##### Theorem 2.

Let fk=ρkd, k{B,C} and f=max(fB,fC). Then we have the following expansion:

(9)σω=k0[I+fBMB0+fCMC0(ω)]+o(f),
where
(10)MC0(ω)ep·eq=C0νq(y)(1β0(ω)k0I+L#,C0)1[νp](y)ds(y),
and
(11)MB0ep·eq=B0(λIK#,B0)1[νp](y)yqds(y)
with
λ=k1+k02(k1k0).

We begin be reviewing properties of periodic layer potentials. Let us define the periodic Green’s function

G#(x)=nZd{0}e2iπn·x4π2|n|2.
Thanks to Poisson’s summation formula, in the sense of distribution, G# satisfies
(12)ΔG#(x)=nZdδ(xn)1.
We write G#(x,y):=G#(xy). Let us introduce the periodic single layer potential, for a Lipschitz domain DY:
S#,D:H1/2(D)Hloc1(RdD)φxDG#(x,y)φ(y)ds(y),
the periodic double layer potential
D#,D:H1/2(D)Hloc1(RdD)φxDG#ν(y)(x,y)φ(y)ds(y),
and the periodic Neumann–Poincaré operator
K#,D:H1/2(D)H1/2(D)φxDG#ν(y)(x,y)φ(y)ds(y),
K#,D:H1/2(D)H1/2(D)φxDG#ν(x)(x,y)φ(y)ds(y).
We review the jump properties of the layer potentials [3].

##### Lemma 3.

We have the following jump relations along the boundary D:

S#,D[φ](x)|+=S#,D[φ](x)|,νS#,D[φ](x)|±=(±12I+K#,D)[φ](x),D#,D[φ](x)|±=(12I+K#,D)[φ](x),νD#,D[φ](x)|+=νD#,D[φ](x)|.
where the subscript ± means fD(x)|±=limt0+fD(x±tν(x)) for xD.

We denote by L#,D the operator φνD#,D[φ]. We write νp=ν·ep on B and C. Using these jump relations, we have the following representation theorem for vp, p{1,,d}.

##### Theorem 4.

We have the following representation for vp:

(13)vp=Cp+S#,B[φ1,p]D#,C[φ2,p],
where Cp is a constant and (φ1,φ2) satisfies the following system:
(14)(λIK#,B)[φ1,p]+νD#,C[φ2,p]=νpon B,(1βk0I+L#,C)[φ2,p]νS#,B[φ1,p]|+=νpon C.

##### Lemma 5.

For any (F,G)H1/2(B)×H1/2(C), the system

(λIK#,B)[φ1]+νD#,C[φ2]=Fon B,(1βk0I+L#,C)[φ2]νS#,B[φ1]|+=Gon C,

##### Proof.

As shown in the Appendix, 1βI+L#,C and λIK#,B are invertible for λ(1/2,1/2]. Moreover, since

νD#,C:H1/2(C)H1/2(B)
and
νS#,B:H1/2(B)H1/2(C)
are compact, the operator
H1/2(Ω)×H1/2(Ω)H1/2(Ω)×H1/2(Ω)(φ1,φ2)((λIK#,B)[φ1]νD#,C[φ2],(1βk0I+L#,C)[φ2]νS#,B[φ1]|+)
is a Fredholm operator. It is therefore sufficient to show that it is injective. Let (φ1,φ2) be such that
(λIK#,B)[φ1]+νD#,C[φ2]=0on B,(1βk0I+L#,C)[φ2]νS#,B[φ1]=0on C.
Let v=S#,B[φ1]D#,C[φ2]. Then v is 1-periodic in every direction, and v is a solution by construction to the following problem:
(15)·(k0(vp(y)+y))=0for yY(BC),·(k1(vp(y)+y))=0for yB,k0ν(vp+(y)+y)=k0ν(vp(y)+y)for yC,vp+(y)vp(y)=β(ω)k0ν(vp+(y)+y)for yC,vp+(y)vp(y)=0for yB,k0ν(vp+(y)+y)=k1ν(vp(y)+y)for yB.
By the uniqueness of the solution to (15) up to a constant, v(x)=c,xY. Then, we have φ1=0 on C and φ2=0 on B because they are equal to the jumps of v (resp. vν) across B (resp. C). This concludes the proof. □

We can now proceed to prove Theorem 4.

##### Proof.

Let (φ1,φ2) be a solution of (14), and let

vp=S#,B[φ1]D#,C[φ2].
Then using the jump relations of the layer potentials, we have that vp is a solution of (8), except that we have not necessarily Yvp=0. We just have to adjust Cp accordingly. □

We now proceed to compute the representation of the effective conductivity.

##### Theorem 6.

We have the following representation for σω:

σω=k0(I+M),
where M=(Mpq)p,q=1d is defined by
(M)pq=Bxpφ1,qdsCνpφ2,qds,p,q{1,,d}.

##### Proof.

We recall the expression of σω in (5):

σωep·eq=k0(δp,q+Yvpν(y)yqds(y)).
Using representation (13), we obtain
Yvpν(y)yqds(y)=YS#,B[φ1,p]ν(y)yqds(y)YD#,C[φ2,p]ν(y)yqds(y)
and
YS#,B[φ1,p]ν(y)yqds(y)=BS#,B[φ1,p]ν|+(y)yqds(y)BS#,B[φ1,p]ν|(y)yqds(y)=Byqφ1,p(y)ds(y).
The same reasoning applies to the second part of the equation:
YD#,C[φ2,p]ν(y)yqds(y)=CD#,C[φ2,p]|+(y)νq(y)ds(y)CD#,C[φ2,p]|(y)νq(y)ds(y)=Cφ2,p(y)νq(y)ds(y).
Therefore,
σωep·eq=k0(δp,q+Yvpν(y)yqds(y))=k0(δp,q+Byqφ1,p(y)ds(y)Cφ2,pνq(y)(y)ds(y)).

We turn to the proof of Theorem 2. We first review asymptotic properties of the periodic Green’s function G#. The following result from [3, Chapter 2] holds.

##### Lemma 7.

We have the following expansion for G#:

G#(x)=G(x)+Rd(x),
where G is the Green function and Rd is a smooth function on Rd and its Taylor expansion at 0 is given by
(16)Rd(x)=Rd(0)12d|x|2+O(|x|4).

Using this expansion, we obtain by exactly the same arguments as those in [3, Chapter 8] the following expansion, which is uniform in zB0,

(λIKB0)[ψB,p](z)=νB0,p(z)+o(1)(1β0k0I+LC0)[ψC,p](z)=νC0,p(z)+o(1),
where KB0 is the standard Neumann–Poincaré operator and LC0 denotes the operator νDC0 associated with the standard double layer potential DC0:
KB0[ϕ](x):=B0Gν(x)(x,y)ϕ(y)ds(y),LC0[ϕ](x):=νC0Gν(y)(x,y)ϕ(y)ds(y).
Therefore, we arrive at the result stated in Theorem 2.

### 3.2.Spectral measure of the tissue

Expansion (9) yields

σω=k0[I+ρBdMB0+ρCdMC0(ω)]+O(ρd)
with
MC0(ω)ep·eq=C0νq(y)(1β0(ω)k0I+LC0)1[νp](y)ds(y).
In order to use the spectral theorem in a Hilbert space, we have to modify the expression of MC0. Let LC01 be the inverse of LC0:H01/2(C0)H01/2(C0). Then we write
(1β0(ω)k0I+LC0)1[νp]=(1β0(ω)k0LC01I+IH1/2)1LC01[νp].
The following result holds.

##### Lemma 8.

LC01I can be extended to a self-adjoint operator L:L2(C0)L2(C0), whose image is a subset of H1/2(C0).

##### Proof.

Let J1:L2(C0)H1/2(C0) and J2:H1/2(C0)L2(C0). Let L=J2LC01J1. Then obviously L extends LC01I and its image is a subset of H1/2(C0). Let us show that it is self-adjoint. Let (φ,ψ)L2(C0)×L2(C0). Let ,L2 and ,H1/2,H1/2 respectively denote the L2-scalar product and the duality pairing between H1/2(C0) and H1/2(C0). We have

L[φ],ψL2=LC01[φ],ψL2=LC01[φ],ψH1/2,H1/2=LC01[ψ],φH1/2,H1/2=LC01[ψ],φL2=L[ψ],φL2,
since LC0 is self-adjoint from H1/2(C0) onto H1/2(C0). □

From this result, we can now proceed. From the spectral theorem, there exists a spectral measure E such that for any zCΛ(L) and for any (φ,ψ)(L2(C0))2,

(17)(Lz+I)1[φ],ψL2=Λ(L)1xz+1φ(x)ψ(x)dE(x),
where Λ(L) denotes the spectrum of L. Let
(18)Fp,q(z)=δp,q+ρBdMB0ep·eq+ρCdΛ(L)1xz+1LC01[νp](x)·νq(x)dE(x),
where δp,q=1 if p=q and δp,q=0 if pq. Therefore, we have
σωep·eqk0[Fp,q(β0(ω)k0)].

Since

limz0F(z)=I+ρBdMB0,
there is no singularity of F in 0. Since 0Λ(L), (17) is valid on a neighborhood of 0.

##### Proposition 9.

Let F=(Fp,q)p,q=1d be defined by (18). Then the following expansion of F in a neighborhood of 0 holds:

(19)Fp,q(z)=k=0ak,p,qzk,
where
a0,p,q=I+ρBdMB0ep·eq,
and
a1,p,q=ρCdνp·νq.

##### Proof.

Identity (19) holds using the analyticity of F in a neighborhood of 0. We also have

a0,p,q=limz0Fp,q(z)=δp,q+ρBdMB0ep·eq.
In order to obtain the next coefficients, we begin by establishing the following limit:
limz0(L+zI)1[νp]=LC0[νp],p=1,2.
Indeed, let φ(z)=(L+zI)1[νp]. Then
φ(z)=1z(νpLφp).
Since the range of L is a subset of H1/2(C0), φ(z)H1/2(C0). Therefore,
φ(z)=LC0[νp]zLC0[φ](z)z0LC0[νp].
This yields
limz01z(Fp,q(z)Fp,q(0))=ρCdνp·νq.

In the following, we write

F(z)=(Fp,q(z))p,q{1,,d},zCΛ(LC0),
and
(20)Ak=(ak,p,q)p,q{1,,d},kN.
Since Fp,q is analytic on CΛ(LC0), the values of ak can be recovered from the values of Fp,q on a subset of C with a limiting point. Therefore, we can reconstruct the values ak,p,q from the measurements of the effective conductivity σω in a band of frequencies ω(ω1,ω2). Further details on this will be provided in the following section.

## 4.Inverse homogenization

### 4.1.Imaging of the anisotropy ratio

The anisotropy ratio (the ratio between the largest and the lowest eigenvalue of the effective conductivity tensor) depends on the frequency [2]. Furthermore, in the general case, the anisotropy orientation (the direction of the effective conductivity tensor eigenvectors) depends also on the frequency. However, in the special case where we have an axis of symmetry of a single inclusion or a cell, the anisotropy orientation is independent of the frequency.

We denote by Od(R):={RRd×d|RTR=1,det(R)=1} the set of rotational matrices. Here, the superscript T denotes the transpose. For convenience, we write R(x):=Rx for xY and R(D):={Rx:xD}. We will need the following covariance result :

##### Lemma 10.

Let ROd(R) and fL2(C0). Then

LC0[fR]R=LC0[f].

##### Proof.

We have, for any xC0,

LC0[fR](R(x))=limh0DC0[fR](R(x)+hν(R(x)))·ν(R(x)).
Moreover,
DC0[fR](R(x))=C0G(R(x)y)·ν(y)f(R(y))ds(y)=C0G(R(x)R(y))·ν(R(y))f(y)ds(y).
Since G is isotropically symmetric,G(R(xy))=R(G(xy)), therefore for any x,yC0,
G(R(x)R(y))·ν(R(y))=R(G(xy))·R(ν(y))=G(xy)·ν(y)
so that
DC0[fR](R(x))=DC0[f](x),xC0.
This in turn implies that
LC0[fR](R(x))=limh0DC0[fR](R(x)+hν(R(x)))·ν(R(x))=limh0DC0[f](x+ν(x))·ν(x)=LC0[f](x).

The following corollary holds immediately.

##### Corollary 11.

Let ROd(R). Then,

MR(C0)=RMC0RT.

Let us begin with the two-dimensional case.

##### Proposition 12.

Let d=2, and (ϵ1,ϵ2) be an orthonormal basis of R2; see Fig. 3. Let ξ be the orthogonal symmetry of axis ϵ1. If ξ(C0)=C0, then

F(z)ϵ1·ϵ2=0,zCΛ(L).

##### Proof.

We have

F(z)ϵ1·ϵ2=ρCdC0(Lz+I)1[ν·ϵ1](x)ν(x)·ϵ2ds(x)=ρCdC0(Lz+I)1[ν·ϵ1](ξ(x))ν(ξ(x))·ϵ2ds(x)=ρCdC0(Lz+I)1[ν·ϵ1](x)ν(x)·ϵ2ds(x)
because ν(ξ(x))·ϵ1=ν(x)·ϵ1 and ν(ξ(x))·ϵ2=ν(x)·ϵ2. Therefore,
F(z)ϵ1·ϵ2=0,zCΛ(L).

We have a similar result in three dimensions. The following proposition holds.

##### Proposition 13.

Let d=3, and (ϵ1,ϵ2,ϵ3) be an orthonormal basis of R3. Let ξ1 (resp. ξ2) be the orthogonal symmetry of axis ϵ1 (resp. ϵ2). If ξ1(C0)=ξ2(C0)=C0, then

F(z)ϵj·ϵk=0,zC,kj{1,2,3}.

##### Proof.

The proof is exactly the same as in the d=2 case and is therefore omitted. □

##### Fig. 3.

A domain presenting a symmetry. In this case, the anisotropy direction is frequency independent.

##### Remark 14.

It is also true that the symmetry axes of B0 correspond to the eigenvectors of the polarization tensor MB0. Therefore, the anisotropy direction of the frequency-independent background can also be recovered as the principal directions of MB0.

##### Remark 15.

Even if each of inclusion and cell has an axis of symmetry, the direction of eigenvectors of the effective conductivity tensor can be frequency dependent. The following numerical test is conducted to show an example of frequency dependency. There are an ellipsoidal inclusion with major axis e1 and minor axis e2 and an ellipsoidal cell with major axis e2 and minor axis e1 in the unit square as shown in Fig. 4(a). For the square domain Y=(12,12)2, each axis length of cell and inclusion is 1/8, and 1/24. The center of ellipsoidal cell and inclusion are (1/3,1/6) and (0,1/3) respectively. The ratio between membrane thickness and size of a cell is 5×103. The conductivity value of medium, membrane, inclusion are 0.5 S/m, 105 S/m, and 1012 S/m respectively. We use (5) to compute the effective conductivity tensor. For the numerical computation, we take advantage of using uj satisfying ·(σuj)=0 in Ω with boundary condition uj(y)|Ω=yj|Ω for y=(y1,y2). Then, vj can be replaced with vj=ujyj. Hence, the eigenvectors of the effective conductivity can be computed and the main direction of anisotropy changes in terms of the frequency as shown in Fig. 4(b).

##### Fig. 4.

(a) shows voltage map with current flows for each y1- and y2-direction current at 104 and 109 Hz. (b) shows eigenvectors of the effective conductivity. Blue arrows represent eigenvectors at frequency ω/2π=104 Hz while red arrows are representing eigenvectors at frequency ω/2π=109 Hz.

### 4.2.Implementation of the inverse homogenization

Following [2], we use the following values:

• The size of cells: 50 μm;

• Ratio between membrane thickness and size of a cell: 0.7×103;

• Medium conductivity: 0.5 S/m;

• Membrane conductivity: 108 S/m;

• Background inclusion conductivity: 107 S/m;

• Membrane permittivity: 3.5×8.85×1012 F/m;

• Frequency band: ω/2π[104;109] Hz.

##### Fig. 5.

Values of β(ω) for ω/2π[104;109].

In this case, we have values of β(ω) for ω/2π[104;109] in Fig. 5. We consider a sample medium as follows: the cells are elliptic in shape, with axes lengths ρCaC and ρCbC, with acbCπ=1. The background is composed of elliptic inclusions, with axes lengths ρBaB and ρBbB, with aBbBπ=1. Their orientation is given by the angles θC and θB respectively.

At each frequency, in order to compute the true effective conductivity given by (5), we perform a finite element computation using freefem++ [7]. Comparison between the true effective conductivity and the expansion from Theorem 9 can be seen in Figs 6 and 7, in the case θB=0 and θC=0, and ρB=ρC=0.1.

##### Fig. 6.

Real part of the effective conductivity.

##### Fig. 7.

Imaginary part of the effective conductivity.

To recover the moments from the effective conductivity, we approximate as a rational function,

Fp,q(z)p0+p1z++pNzNq0+q1z++qNzN
for some NN. Such an approximation of F is called a Padé approximation of F. Then we approximate the moments by the following values:
a˜0,p,q=p0q0,a˜1,p,q=p1q0q1p0q02.

Numerically, this is done as a simple least square inversion: the coefficients of the polynomials P(z)=p0+p1z++pNzN and Q(z)=q0+q1z++qNzN are computed to minimize the quantity

k=1K|Fp,q(zk)P(zk)Q(zk)|2,
where z1,,zK are the frequency values where F is measured.

We now consider a toy example where C is an ellipse in R2. In this case, if λ1 and λ2 are the eigenvalues of A1 defined by (20) for k=1, the ratio r:=λ2/λ1 is independent of the volume fraction and is given by

(21)r=02πb2cos2(t)b2cos2(t)+a2sin2(t)dt02πa2sin2(t)b2cos2(t)+a2sin2(t)dt=ba02πcos2(t)cos2(t)+a2b2sin2(t)dt02πsin2(t)b2a2cos2(t)+sin2(t)dt.

Since the right-hand side of (21) can be regarded as a function of a/b, the anisotropy ratio a/b can be easily obtained by solving (21) with the known value r. In Fig. 8 (resp. in Fig. 9), we illustrate the reconstruction of the ratio r using the Padé approximation of F as a function of the anisotropy ratio a/b compared to its theoretical value given by the preceding formula in the case where there is no inclusion B (resp. with an inclusion B with ρB=0.1). As we can see, the reconstruction is almost perfect in the case where there is no inclusion, and there is a slight bias induced by the inclusion B.

##### Fig. 8.

Reconstruction of r when there is no inclusion B.

##### Fig. 9.

Reconstruction of r when there is an inclusion B with ρB=0.1.

After recovering the anisotropy ratio a/b, we can recover the volume fraction ρC from the product of λ1,λ2 of the eigenvalues of A1. Indeed, we have

λ1λ2=ρC4ab02πcos2(t)cos2(t)+a2b2sin2(t)dt02πsin2(t)b2a2cos2(t)+sin2(t)dt=ρC4π02πcos2(t)cos2(t)+a2b2sin2(t)dt02πsin2(t)b2a2cos2(t)+sin2(t)dt.
Table 1 presents numerical reconstruction of the volume fraction ρC using the preceding formula, with an anisotropy ratio equal to 2.

##### Table 1

Reconstructed values of ρC with anisotropy ratio of 2

 Values of ρC 0.01 0.02 0.03 0.05 0.1 0.2 0.3 Reconstructed value 0.0098 0.0196 0.0294 0.0491 0.0981 0.1963 0.2945

To reconstruct the angle of the inclusions, we simply use the orientation of the eigenvalues of the moments of A0 for B and A1 for C. This is illustrated by results in Fig. 10 when both B and C are ellipses of anisotropy ratio 2 and with ρB=ρC=0.1.

##### Fig. 10.

Reconstruction of the orientation of the inclusions B and C.

## Appendices

### AppendixSpectrum of some periodic integral operators

Let CRd be a C1+α-domain for some α>0. It is known that the non periodic operator λIKC is invertible on H1/2 for λ(12,12] [4,6]. The positivity of LC [12, Section 3.3] also implies that λI+LC:H1/2H1/2 is invertible for λ>0. We extend these results to the case of periodic Green’s function.

##### Theorem 16.

For any λ>0, the operator λI+L#,C:H1/2(C)H1/2(C) is invertible.

##### Proof.

We first show that the operator L#,C is a Fredholm operator. Note that, L#,C=LC+R where R is an integral operator with a smooth kernel and is therefore compact. Moreover, since LC has a dimension 1 kernel and image, it is a Fredholm operator. Therefore, L#,C is Fredholm. Now we show that L#,C is positive semi-definite, and the result will follow from the Fredholm alternative. Since

L#,C[φ],ψL2=S#,C[curlCφ],curlCψL2
for any φ,ψH1/2(C), we just have to show that S#,C is negative semi-definite. From the expression (12) for G#, we compute, for any φL2(C),
S#,C[φ],φL2=nZd{0}CCe2iπn·(xy)4π2|n|2φ(x)φ(y)ds(x)dS(y)=nZd{0}(Ce2iπn·y2π|n|φ(y)dS(y))(Ce2iπn·x2π|n|φ(x)ds(x))=nZd{0}|Ce2iπn·y2π|n|φ(y)ds(y)|20.
Therefore, S#,C is negative semi-definite, which concludes the proof. □

##### Theorem 17.

For λ(12,12], the operator λIK#,C is invertible on H1/2(C).

##### Proof.

Since λIKC is invertible, K#,CKC is a compact operator [3], λIK#,C is a Fredholm operator and it is enough to show that it is one-to-one. The proof goes exactly as in [4]. Let us assume that λIK#,C is not one-to-one. Then there exists some fH1/2(C) such that

(λIK#,C)[f]=0.
Let us write
(λIK#,C)[f]=(λ12)f+(12IK#,C)[f].
Since (12IK#,C)[f],1L2=0, we have f,1L2=0. Let u=S#,C[f]H1(YC). Let
A=C|u(x)|2dxandB=YC|u(x)|2dx.
Then A0 or B0 since f is not identically zero. Then by Green’s formula together with the jump formulas, we have
A=(12I+K#,C)[f],S#,C[f]L2andB=(12I+K#,C)[f],S#,C[f]L2.
Since (λIK#,C)[f]=0, we have β=12BAB+A. We have therefore a contradiction: we have |β|12 since A,B0. Therefore, β=12 which implies that B=0. Therefore, u is constant in RdnZd{C+n}. Since u is continuous across C, u is harmonic on C and is constant on C, and by uniqueness of the Dirichlet problem on C, u is constant on C. Therefore,
f=νS#,C[f]|+νS#,C[f]|=0,