Mostafa Bendahmane
email: mostafa.bendahmane@u-bordeaux.fr

Elmahdi Erraji
email: erraji0elmahdi@gmail.com

Fahd Karami
email: fa.karami@uca.ma

A 3D REACTION-DIFFUSION SYSTEM DESCRIBING CALCIUM DYNAMICS IN CARDIAC CELL *

Keywords: Mathematics Subject Classification. 92C50, 65M60, 65N12, 92C40 Finite element, Calcium dynamics, Anomalous diffusion, Convergence, Interface problem, Medical modeling

We are interested in modeling the interaction of calcium dynamics in a medium including sarcolemma and sarcoplasmic reticulum. The governing equations consist of a nonlinear reactiondiffusion system representing the various calcium fluxes and theirs buffers in the two media. We address the question of existence of weak solutions by using a fixed-point approach. We propose a finite element method for this system, we establish the existence of the discrete solution, and we show that the discrete solution generated by the given scheme converges to the corresponding weak solution for the model studied. Finally, we give some 2D and 3D numerical examples to our model.

Introduction

The study of calcium dynamics in cardiac cell is a priority to reveal the mysterious causes of cell failure. One of the principal causes of heart failure is the perturbation of calcium induced-calcium release (CICR) cycle of a cardiac cell (see for e.g. [START_REF] Marks | Calcium and the heart : a question of life and death[END_REF]). The calcium cycle in cardiac cell is a complicated process, since it involves different mechanisms such as swapping between the cytoplasm and sarcoplasmic reticulum through different type of channels and interactions with biochemical species in the cell medium. On the cellular scale, the perturbation of CICR could be explained by pathological behavior of ionic channels especially Ryanodine receptor (RyR). The scientific community is convinced, based on strong evidence (we refer to [START_REF] Vest | Defective cardiac ryanodine receptor regulation during atrial fibrillation[END_REF]) that a vast majority of pathological heart cases is due to RyR irregularity. Which induces more questions than answers about understanding the causes and designing an effective therapeutic strategy on the cellular scale. This would require a larger collaboration between various scientific fields including physics, chemistry, biology and mathematics.

In cardiac and other muscle cells, calcium signaling is the main factor to mediate excitation-contraction coupling (ECC). Its cycle lies between the action potential and muscle contraction. The ECC is a consequence of a propagating action potential through muscle cell membranes. At the cellular scale, action potential propagates along the cell membrane including the T-tubules which forms a tubular cavity inside the cell. This action potential triggered voltage opened a channel (VOC) or L-type channels located on the cell membrane surface. As a result, a passive flux of Ca 2+ crosses the L-type due to the large difference of calcium concentration between the extracellular ([Ca 2+ ] e ≈ 1mM ) and the intracellular ([Ca 2+ ] i ≈ 0.1µM ) medium. The rising Ca 2+ concentration triggered the calcium-induced calcium release process (CICR). It begins by triggering a ryanodine receptor located on sarcoplasmic reticulum, this allows the Ca 2+ influx into intracellular medium from sarcoplasmic reticulum, where Ca 2+ is stored at elevated concentration ([Ca 2+ ] sr ≈ 1.3mM , see [START_REF] Bers | Excitation-contraction coupling and cardiac contractile force[END_REF] for this detail), not to mention the calcium concentration maintained by the heavy buffering. This leads to the simultaneous opening of the other ionic channel. Then, the influx of Ca 2+ ion diffuses through myofibrils to cause mechanical activity.

In literature, many authors attempt to build a primal model of calcium dynamics. We mention, for example the works of Peskoff et al. in [START_REF] Langer | Calcium Concentration and Movement in the Diadic Cleft Space of the Cardiac Ventricular Cell[END_REF][START_REF] Peskoff | sarcolemmal calcium binding sites in heart: II. Mathematical model for diffusion of calcium released from the sarcoplasmic reticulum into the diadic region[END_REF], their first work establishes a primal model of calcium distribution evolution in diad neighbor. A more complete calcium model is established by extending the previous model to the whole cell and considering interaction with local buffers in the medium (see [START_REF] Picht | Dynamic calcium movement inside cardiac sarcoplasmic reticulum during release[END_REF] for more details). In passing, we want to mention the work of Sobie et al. in [START_REF] Sobie | Termination of Cardiac Ca2+ Sparks: An Investigation Mathematical Model of Calcium-Induced Calcium Release[END_REF], where the authors studied the dynamics of Ca 2+ under new consideration on RyR Stochastic behavior, period of release and there distribution over space. Note that these studies of spatial calcium dynamics doesn't take in consideration the complexity of cell medium (for example the presence of mitochondria in the cell and the heterogenerous of the media). This will affect the calcium diffusion and produce some serious needs in the modeling part. Recently, some new spatial calcium models are based on a transmission defined by interaction in a local position of ionic channels (see for e.g. [START_REF] Chai | Towards simulation of subcellular calcium dynamics at nanometre resolution[END_REF][START_REF] Langer | Calcium Concentration and Movement in the Diadic Cleft Space of the Cardiac Ventricular Cell[END_REF]). These models can be interpreted as a transmission boundary value problem and can be studied by using domain decomposition techniques. From the analytical and the numerical perspectives, there are many results that has been established in this context (see for e.g. [START_REF] Flauraud | Optimized interface conditions in domain decomposition methods for problems with extreme contrasts in the coefficients[END_REF][START_REF] Lagnese | Dynamic Domain Decomposition in Approximate and Exact Boundary Control in Problems of Transmission for Wave Equations[END_REF][START_REF] Mayer | Interface handling for three-dimensional higher-order XFEM-computations in fluidstructure interaction[END_REF]). Note that the previous spatial calcium models are governed by a linear diffusion. We know that the solution of a linear diffusion equation admits an infinite speed of propagation whereas the solution of the nonlinear diffusion equation shows a finite speed of propagation. Moreover, from the application point of view it is not realistic to have the property of infinite speed of propagation because the calcium can not diffuse infinitely. In this paper, we establish the mathematical model of species dynamics in cardio-myocyte with nonlinear diffusion. Note that strong experimental evidence proves the complexity of calcium diffusion (see for e.g. [START_REF] Donahue | Free diffusion coefficient of ionic calcium in cytoplasm[END_REF]). Theoretically, the ionic nature of Ca 2+ as well as the complexity of the cytosolic and sacroplasmic medium imposes some modifications on the previous models (we refer here to the works [START_REF] Pods | Electrodiffusion Models of Neurons and Extracellular Space Using the Poisson-Nernst-Planck Equations-Numerical Simulation of the Intra and Extracellular Potential for an Axon Model[END_REF][START_REF] Rubinstein | Electro-Diffusion of ions. 1st edition[END_REF]). Moreover, advanced microscopical imaging technique shows that the cytoplasm of mammalian cardiac myocytes seems to be highly heterogeneous and anisotropic due to the dense presence of some organelles. We mention, for example the presence of mitochondria and Golgi complex on the micron-scale. A zoom into a smaller scale, clarify the presence of a network of elements such as microtubule and actin filaments. Hence, based on strong experimental evidence shown in [START_REF] Gal | Experimental evidence of strong anomalous diffusion in living cells[END_REF], the authors highly recommend to consider an anomalous diffusion in the cell medium. Based on the same arguments between others [START_REF] Severs | Remodelling of gap junctions and connexin expression in diseased myocardium[END_REF], a recent work [START_REF] Hurtado | Computational modeling of non-linear diffusion in cardiac electrophysiology: A novel porous-medium approach[END_REF] has considered a porous medium model describing cardiac-electro-physiology.

Our main contribution in this context is to generalize the previous models to a nonlinear model with non constant diffusion coefficients. This paper is devoted to a rigorous study of existence and uniqueness of L ∞ bounded solutions of the proposed model. Moreover, we establish a fully discrete scheme in finite element framework for spatial discretization, whereas the first order backward Euler scheme is used for time discretization. Then, we move forward to prove the convergence result of the proposed scheme. We want to mention here (in the context of finite element method) that we cannot obtain the maximum principle as in the continuous case because the positive and negative cuts test functions may not be admissible. For the convergence of the scheme, we will assume that the initial conditions are only bounded in L 2 . The plan of this paper is as follows: In Section 1 we propose the derivation of the calcium model based on conservation law and interaction between species. Section 2 is devoted to the existence and uniqueness of weak solution (defined in Definition 2.1 below) by using the Schauder fixed point theorem and the passage to the limit that transforms the regularized system into the original problem. We introduce in Section 3 some notations for the finite element method and we present our scheme. The proof of the convergence result (of the scheme) is divided into a priori estimates, compactness for discrete solution and convergence to a weak solution. Finally, we give in Section 4 some numerical examples to our model.

Calcium transmission model

We devote this section to a brief derivation of the calcium model of cardiac tissue. As principal references for this model with constant diffusion we found [START_REF] Chai | Towards simulation of subcellular calcium dynamics at nanometre resolution[END_REF]. Before starting our calcium dynamics modeling, we refer to Figure 1 for a good description of the cytosolic and sarcoplasmic reticulum domains.

. Description of calcium transmission dynamics

In both sub-domains, calcium behaves in a similar manner. To derive the partial differential equation controlling calcium behavior, we consider an arbitrary domain Ω. By the conservation law, we have

d dt Ω c dV = Ω R dV - ∂Ω J • η dA, ( 1
)
where c is the calcium concentration, J its flux trough ∂Ω, R its source inside Ω due to interaction of Ca 2+ with other species and η is the normal vector on ∂Ω. Now to conclude this equation, we use Stokes Ostradski formula

∂Ω J • η dA = Ω div J dV. (2)
Combining ( 1) and ( 2), we arrive

d dt Ω c dV = Ω R dV - Ω div J • η dV. (3)
This implies

∂ t c = R -div J in (0, T ) × Ω. ( 4
) Thanks to Fick law, the calcium moves from high concentration to low concentration in cell medium. So the flux J is directed in the opposite direction as the concentration gradient. Moreover, the diffusion coefficient depend on solvent concentration. In summary, calcium flux can be written as

J = -D c (c)∇c, (5)
where D c is a diffusion-dependent function of c. Hence, we deduce :

∂ t c = div (D c (c)∇c) + R c in (0, T ) × Ω, (6)
and similarly the equation of Ca 2+ in SR is:

∂ t s = div (D s (s)∇s) + R s in (0, T ) × Ω, (7)
where D s is a diffusion-dependent function of s. We will define the sources R c and R s explicitly in the following sub-sections.

Description of buffers dynamics

In this subsection, we establish the equations governing by the buffers (recall that the buffers are proteins that interact with calcium to produce another molecule). We assume that the total concentration of the free buffer [buffer] and the new buffer [buffer-Ca 2+ ] (resulting from the interaction of the free buffer with Ca 2+ ) is a constant in time. The chemical equation describing this phenomenon is:

Buffer + Ca 2+ kon k of f Buffer-Ca 2+ , (8)
where k of f and k on are two rates coefficients.

Next we let b the concentration of buffer-Ca 2+ , B tot the total concentration of the buffer and B the concentration of buffer only . An application of the law of mass action (see for e.g [START_REF] Kenakin | The mass action equation in pharmacology[END_REF]), we obtain

∂ t b = R(c, b) = k on cB -k of f b. (9)
We notice that B = B tot -b so the equation becomes:

∂ t b = R(c, b) = k on c(B tot -b) -k of f b. (10)
Each buffer (for more information, check [START_REF] Perin | An Essential Guide to Cardiac Cell Therapy[END_REF]) is described by a conservation law similar to [START_REF] Anaya | On a vorticity-based formulation for reaction-diffusion-Brinkman systems[END_REF]. So we conclude a similar PDE. The buffer list living in the cytosol is :

• ATP or nucleoside triphosphate used in the cells as a coenzyme often called the molecular currency of intracellular energy transfer.

∂ t b 1 = D b1 ∆b 1 + R 1 (c, b 1 ), (11)
where D b1 is a diffusion coefficient and the source term

R 1 (c, b 1 ) := k on 1 c(B 1 -b 1 ) -k of f b 1 . • CMDN or calmodulin is a multi-functional intermediate calcium-binding messenger protein. ∂ t b 2 = D b2 ∆b 2 + R 2 (c, b 2 ), ( 12
)
where D b2 is a diffusion coefficient and the source term

R 2 (c, b 2 ) := k on 2 c(B 2 -b 2 ) -k of f b 2 .
• Fluo-4 is used in the life sciences generally as a non destructive way of tracking or analyzing biological molecules such as Ca 2+ .

∂ t b 3 = D b3 ∆b 3 + R 3 (c, b 3 ), (13)
where D b3 is a diffusion coefficient and the source term

R 3 (c, b 3 ) = k on 3 c(B 3 -b 3 ) -k of f b 3 .
• TRPN or troponin is a protein of muscle that together with tropomyosin forms a regulatory protein complex controlling the interaction of actin and myosin and when combined with calcium ions permits muscular contraction. Similarly to Fluo-4, the protein TRPN satisfies the following equation

∂ t b 4 = D b4 ∆b 4 + R 4 (c, b 4 ), (14)
where D b4 is a diffusion coefficient and the source term R 4 (c, b 4 ) = k on 4 c(B 4 -b 4 ) -k of f b 4 . We consider only one buffer living in SR :

• CSQN or calsequestrin is a calcium binding protein. It helps hold calcium in the citerna of the sarcoplasmic reticulum after a muscle contraction. Moreover, the protein CSQN is governed by the equation

∂ t b 5 = D b5 ∆b 5 + R 5 (s, b 5 ), ( 15
)
where D b5 is a diffusion coefficient and the source term R 5 (s, b 5 ) = k on 5 s(B 5 -b 5 ) -k of f b 5 . The interaction between buffers and calcium define the source terms R c and R s mentioned in ( 6) and [START_REF] Chen | Calcium waves initiating from the anomalous subdiffusive calcium sparks[END_REF].

R c = - 4 i=1 R i , R s = -R 5 .

Boundary condition

In cardiac cell's, the calcium concentration in SR is considerably high. The activation of ionic channels and Ca 2+ low concentration in cytosol induces a calcium flow throw RyR proportional to its difference between the two sub-domains: where η c (Resp. η s ) is the outward normal vector with respect to Ω cyt (Resp. Ω sr ), Γ ryr is subset of ∂Ω sr ∪∂Ω cyt and Γ ryr,T := Γ ryr × (0, T ). Since the SR and cell membrane are impermeable. Then, there is no flux through it. This can be implemented as an homogeneous Neumann boundary condition D s (s)∇s • η s = 0 on Γ cyt iso,T := Γ cyt iso × (0, T ), and

J • η = ±g(s -
D c (c)∇c • η c = 0 on Γ sr iso,T := Γ sr iso × (0, T ),
where Γ i iso := ∂Ω i -Γ ryr for i = cyt, sr. The flux of buffers is null through cell membrane

∇b i • η c = 0, on Γ T := Γ × (0, T ), for i = 1, 2, 3, 4,
and

∇b 5 • η s = 0 on Γ sr T := Γ sr × (0, T ).
Finally, combining the previous equations, we get the following nonlinear reaction-diffusion system

                                               ∂ t c = div (D c (c)∇c) - 4 i=1 R i (c, b i ) in Ω cyt,T := Ω cyt × (0, T ), ∂ t b i = D bi ∆b i + R i (c, b i ) in Ω cyt,T for i = 1, . . . , 4, ∂ t s = div (D s (s)∇s) -R s (s, b 5 )
in Ω sr,T := Ω sr × (0, T ),

∂ t b 5 = D b5 ∆b 5 + R 5 (s, b 5 ) in Ω sr,T , D c (c)∇c • η c = 0 on Γ cyt iso,T := Γ cyt iso × (0, T ), D s (s)∇s • η s = 0 on Γ sr iso,T := Γ sr iso × (0, T ), ∇b i • η c = 0 on Γ T := Γ × (0, T ), i = 1, . . . , 4, ∇b 5 • η c = 0 on Γ sr T := (Γ i ∪ Γ ryr ) × (0, T ), D s (s)∇s • η s = -g(s -c) on Γ ryr,T := Γ ryr × (0, T ), D c (c)∇c • η c = g(s -c) on Γ ryr,T := Γ ryr × (0, T ), (c, b i )(•, 0) = (c 0 , b i,0 )(•) in Ω cyt and (s, b 5 )(•, 0) = (s 0 , b 5,0 )(•) in Ω sr . (16)
In this system, the nonlinear diffusivities D c , D s are assumed continuous, positive, coercive, and there exist constants 0 < D min < D max :

D s , D c continuous, 0 < D min ≤ D s (m), D c (m) ≤ D max < ∞ for m ∈ R, and 0 < D min ≤ D bi where D bi is a constant diffusion for i = 1, . . . , 4. (17)
Remark 1. Because of the membrane impermeability, we considered in our model the Neumann boundary condition on buffers equations. For calcium dynamics in both sub-domain, the cell membrane contains the ionic channel RyR that allows only the transmission of Ca 2+ . Hence, we have restricted the transmission boundary condition only for calcium equations on Γ ryr . Moreover, we want to mention that the initial data b i,0 of the buffer should be less than its total concentration B i (Total concentration with or without calcium) for i = 1, . . . , 5. In the paper [START_REF] Centor | Diffusion of troponin testing in unstable angina patients: adoption prior to guideline release[END_REF], the author has studied the diffusivity of some buffers (In our case b 4 and b 5 ). Without loss of generality, in this work we consider the small diffusion coefficient in troponin (b 4 ) and calsequestrin (b 5 ) equations.

2.

Well-posedness (existence and uniqueness) of the weak solution: Continuous case.

In this section, we introduce the variational formulation of the problem [START_REF] Lagnese | Dynamic Domain Decomposition in Approximate and Exact Boundary Control in Problems of Transmission for Wave Equations[END_REF]. Then, we look after existence results.

Definition 2.1 (L 2 bound solution). A weak L 2 bound solution of (16) is a seven component (c, b 1 , b 2 , b 3 , b 4 , s, b 5 ) such that c, b i ∈ L ∞ (0, T ; L 2 (Ω cyt )) ∩ L 2 (0, T ; H 1 (Ω cyt )) and s, b 5 ∈ L ∞ (0, T ; L 2 (Ω sr )) ∩ L 2 (0, T ; H 1 (Ω sr )), THE EDPSMATH DOCUMENTCLASS USERS GUIDE VERSION 2 ∂ t c, ∂ t b i ∈ L 2 (0, T ; H 1 (Ω cyt ) ), ∂ t s, ∂ t b 5 ∈ L 2 (0, T ; H 1 (Ω sr ) )
, and satisfying Our first main result is the following existence theorem for weak solutions. To prove theorem 2.3 we first prove existence of solutions to the approximate problem [START_REF] Louch | Control of Ca2+ Release by Action Potential Configuration in Normal and Failing Murine Cardiomyocytes[END_REF] below by applying the Schauder fixed-point theorem (in an appropriate functional setting). Then we send the regularization parameter to zero to produce a weak solution of the original system (16) as the limit of a sequence of such approximate solutions. Convergence is achieved by means of a priori estimates and compactness arguments.

T 0 ∂ t c, ϕ c cyt dt + Ω cyt,T D c (c)∇c • ∇ϕ c dx dt - Γ ryr,T g(s -c)ϕ c dσdt = - Ω cyt,T 4 i=1 R i (c, b i )ϕ c dx dt, T 0 ∂ t s, ϕ s sr dt + Ω sr,T D s (s)∇s • ∇ϕ s dx dt + Γ ryr,T g(s -c)ϕ s dσdt = - Ω sr,T R 5 (s, b 5 )ϕ s dx dt, T 0 ∂ t b i , ϕ i cyt dt + D bi Ω cyt,T ∇b i • ∇ϕ i dx dt = Ω cyt,T R i (c, b i )ϕ i dx dt, T 0 ∂ t b 5 , ϕ 5 sr dt + D b5 Ω cyt,T ∇b 5 • ∇ϕ 5 dx dt = Ω sr,T R 5 (c, b 5 )ϕ 5 dx dt, (18)
for all ϕ c , ϕ i ∈ L 2 (0, T ; H 1 (Ω cyt )), ϕ s , ϕ 5 ∈ L 2 (0, T ; H 1 (Ω sr )) for i = 1, . . .

Existence result to the regularized problem: the fixed-point method

In this subsection we prove, for each fixed ε > 0, the existence of solutions to the approximate problem [START_REF] Louch | Control of Ca2+ Release by Action Potential Configuration in Normal and Failing Murine Cardiomyocytes[END_REF] below, by applying the Schauder fixed-point theorem. For this purpose, we introduce the following closed subset of the Banach space

A = {(c, s) ∈ H : 0 ≤ c(t, xc), s(t, xs) ≤ R for a.e. (t, xc) ∈ Ω cyt,T and (t, xs) ∈ Ω sr,T },
for some constant R > 0. Herein, H := L 2 (Ω cyt,T ) × L 2 (Ω sr,T ). For a given nonnegative function Û = (ĉ, ŝ) such that ĉ ∈ L 2 (Ω cyt,T ) and ŝ ∈ L 2 (Ω sr,T ), we solve the approximate system

                                               ∂ t c = div (D c (ĉ)∇c) - 4 i=1 R i,ε (c + , b + i ) in Ω cyt,T , ∂ t b i = D bi ∆b i + R i,ε (ĉ + , b + i ) in Ω cyt,T , for i = 1, . . . , 4, ∂ t s = div (D s (ŝ)∇s) -R s,ε (s + , b + 5 ) in Ω sr,T , ∂ t b 5 = D b5 ∆b 5 + R 5,ε (ŝ + , b + 5 ) in Ω sr,T , D c (c)∇c • η c = 0 on Γ cyt iso,T , D s (s)∇s • η s = 0 on Γ sr iso,T , ∇b i • η c = 0 on Γ T , for i = 1, . . . , 4 ∇b 5 • η s = 0 on Γ sr T , D s (s)∇s • η s = -g(s -c) on Γ ryr,T , D c (c)∇c • η c = g(s -c) on Γ ryr,T , (c, b i )(•, 0) = (c 0 , b i,0 )(•) in Ω cyt and (s, b 5 )(•, 0) = (s 0 , b 5,0 )(•) in Ω sr , (19)
defined for each ε > 0. Herein,

R κ,ε = R κ 1 + ε |R κ | for κ = 1, . . . , 5. (20)
Next, we introduce the map F defined by

F : A -→ A Û -→ U,
where U = (c, s) solves the approximate system [START_REF] Louch | Control of Ca2+ Release by Action Potential Configuration in Normal and Failing Murine Cardiomyocytes[END_REF]. The following lemma shows some estimates that will guarantee the continuity and compactness of F below.

||c|| L ∞ (0,T ;L 2 (Ωcyt)) + ||c|| L 2 (0,T ;H 1 (Ωcyt)) ≤ C, ||b i || L ∞ (0,T ;L 2 (Ωcyt)) + ||b i || L 2 (0,T ;H 1 (Ωcyt)) ≤ C for i = 1, . . . , 4, ||s|| L ∞ (0,T ;L 2 (Ωsr)) + ||s|| L 2 (0,T ;H 1 (Ωsr)) ≤ C, ||b 5 || L ∞ (0,T ;L 2 (Ωsr)) + ||b 5 || L 2 (0,T ;H 1 (Ωsr)) ≤ C. (21)
Proof. For the proof it suffices to use assumptions on the initial conditions and the definition of functions R i (for i = 1, . . . , 5) combined with the general result for parabolic equations (see for e.g. for more details [START_REF] Ladyzhenskaya | Linear and quasi-linear equations of parabolic type[END_REF]).

To prove the existence and uniqueness of ( 16), we apply now the Schauder fixed point theorem to the map F.

The goal is to prove that such map has a fixed-point. First, we notice that F is well defined. To satisfy the hypothesis of the fixed point theorem, we show that F is continuous and compact mapping. Firstly, we prove the continuity of F. Let (ĉ n , ŝn ) be a sequence in A and (ĉ, ŝ) ∈ A be such that

(ĉ n , ŝn ) → (ĉ, ŝ) in L 2 (Ω cyt,T ) × L 2 (Ω sr,T ), as n → ∞. Let us then define (c n , s n ) = F(ĉ n , ŝn ),
i.e., (c n , s n ) is the solution of ( 19) associated with the sequence (ĉ n , ŝn ). The objective is to show that:

(c n , s n ) converges to F(ĉ, ŝ), in L 2 (Ω cyt,T ) × L 2 (Ω sr,T
). We start with the following lemma:

Lemma 2.5. Let (c n , b 1,n , b 2,n , b 3,n , b 4,n , s n , b 5,n
) n be the solution to problem [START_REF] Louch | Control of Ca2+ Release by Action Potential Configuration in Normal and Failing Murine Cardiomyocytes[END_REF]. Then (i) There exists a constant C > 0 such that for i = 1, . . . , 4

0 ≤ c n (x c , t), b i,n (x c , t), s n (x s , t), b 5,n (x s , t) ≤ C, for a.e. (x c , t) ∈ Ω cyt,T and (x s , t) ∈ Ω sr,T . (ii) The sequence (c n , b 1,n , b 2,n , b 3,n , b 4,n , b 5,n , s n ) n is bounded in L 2 (0, T ; H 1 (Ωcyt)) ∩ L ∞ (0, T ; L 2 (Ωcyt)) 5 × L 2 (0, T ; H 1 (Ωsr)) ∩ L ∞ (0, T ; L 2 (Ωsr)) 2 . (iii) The sequence (c n , b 1,n , b 2,n , b 3,n , , b 4,n , b 5,n s n ) n is relatively compact in (L 2 (Ω cyt,T )) 5 × (L 2 (Ω sr,T )) 2 .
Proof. To prove (i), we multiply the first and the third equations in ( 19) by -c - n and -s - n (where α -:= (α -|α|)/2) and we integrate over Ω cyt and Ω sr , respectively. The result is 1 2

d dt Ωcyt c - n 2 dx - Ωcyt D c (ĉ n )∇c n • ∇c - n dx + Γryr g(s n -c n )c - n dσ dt = 4 i=1 Ωcyt R i,ε (c + n , b + i,n )c - n dx ≤ 0, and
1 2 d dt Ωsr s - n 2 dx - Ωsr D s (ŝ n )∇s n • ∇s - n dx - Γryr g(s n -c n )s - n dσ dt = Ωsr R s,ε (c + n , b + i,n )s - n dx ≤ 0.
Using this and (17), we deduce 1 2

d dt ||c - n || 2 L 2 (Ωcyt) + 1 2 d dt ||s - n || 2 L 2 (Ωsr) + Γryr g(s -c)(c - n -s - n )dσ dt ≤ 0.
Observe that the function (•) -is a decreasing function. This implies

Γ ryr,T g(s n -c n )(c - n -s - n )dσ dt ≥ 0.
Therefore we obtain 1 2

d dt ||c - n || 2 L 2 (Ωcyt) + 1 2 d dt ||s - n || 2 L 2 (Ωsr) ≤ 0. Since the initial condition (c 0 , s 0 ) is nonnegative we conclude c n (t, x c ), s n (t, x s ) ≥ 0 for a.e. (t, x c ) ∈ Ω cyt,T and (t, x s ) ∈ Ω sr,T . (22)
Reasoning similarly we get the nonnegativity of b i,n for i = 1, . . . , 5.

In the next step we prove first (recall that B i is the total concentration of the buffer)

b i,n (t, x c ) ≤ B i and b 5,n (t, x s ) ≤ B 5 for a.e. (t, x c ) ∈ Ω cyt,T and (t, x s ) ∈ Ω sr,T , (23)
for i = 1, . . . , 4. For this, we multiply (b i,n -B i ) + (where α + := (α + |α|)/2) by the corresponding equation of b i,n in [START_REF] Louch | Control of Ca2+ Release by Action Potential Configuration in Normal and Failing Murine Cardiomyocytes[END_REF] and we integrate over Ω cyt for i = 1, . . . , 4 and Ω sr for i = 5, respectively, to deduce from the definition of R i,ε

1 2 4 i=1 d dt ||((b i,n -B i ) + || 2 L 2 (Ω cyt ) + 1 2 d dt ||((b 5;n -B 5 ) + || 2 L 2 (Ωsr ) + 3 i=1 ||∇(b i,n -B i ) + || 2 L 2 (Ω sr,T ) = 4 i=1 Ω cyt,T R i,ε (ĉn, b i,n )(b i,n -B i ) + dx dt + Ωsr R 5,ε (ŝn, b 5,n )(b 5,n -B 5 ) + dx dt ≤ 0.
Using this and the Gronwall lemma we get [START_REF] Perin | An Essential Guide to Cardiac Cell Therapy[END_REF]. To complete the prove, let us define a super solution to the equations of calcium :

   d dt M (t) = - 5 i=1 R i (M (t), 0) = - 5 i=1 k on i M (t)B i , M (0) = ||c 0 || L ∞ (Ωcyt) + ||s 0 || L ∞ (Ωsr) . (24)
Now, we multiply the first and the third equations in ( 19) by (c n -M (t)) + and (s n -M (t)) + , respectively, and we integrate over Ω cyt and Ω sr , respectively, and we sum the resulting equations to obtain

1 2 d dt ||(cn -M (t)) + || 2 L 2 (Ω cyt ) + 1 2 d dt ||(sn -M (t)) + || 2 L 2 (Ωsr ) + ||∇(cn -M (t)) + || 2 L 2 (Ω cyt,T ) + ||∇(sn -M (t)) + || 2 L 2 (Ω sr,T ) + Γryr g(sn -cn)((sn -M (t)) + -(cn -M (t)) + )dσ dt = 4 i=1 Ω cyt,T (R i,ε (M (t), 0) -R i,ε (cn, b i,n )) (cn -M (t)) + dx dt + Ωsr (Rs,ε(M (t), 0) -Rs,ε(sn, b 5,n )) (sn -M (t)) + dx dt = 4 i=1 Ω cyt,T (R i,ε (M (t), 0) -R i,ε (cn, 0) + R i,ε (cn, 0) -R i,ε (cn, b i,n )) (cn -M (t)) + dx dt + Ωsr (Rs,ε(M (t), 0) -Rs,ε(sn, 0) + Rs,ε(sn, 0) -Rs,ε(sn, b 5,n )) (sn -M (t)) + dx dt. (25)
By monotonicity of (• -M (t)) + we have Γ ryr,T g(s n -c n )((s n -M (t)) + -(c n -M (t)) + )dσ dt ≥ 0. For the right hand terms, we assume that c n ≥ M (t) and s n ≥ M (t). Then, by monotonicity of R i,ε we conclude that :

1 2 d dt ||(c n -M (t)) + || 2 L 2 (Ωcyt) + 1 2 d dt ||(s n -M (t)) + || 2 L 2 (Ωsr) ≤ 0. (26)
Therefore, we have c n (t, x c ), s n (t, x s ) ≤ M (t) for all (t, x c ) ∈ Ω cyt,T and (t, x s ) ∈ Ω sr,T . Moreover, an application of the Gronwall lemma on ( 24) we obtain

M (t) ≤ M (0)e -tC ≤ ||c 0 || L ∞ (Ωcyt) + ||s 0 || L ∞ (Ωsr) with C = 5 i=1 B i k on i .
To prove (ii), we multiply the first, the second, the third, the fourth and the fifth equation in [START_REF] Louch | Control of Ca2+ Release by Action Potential Configuration in Normal and Failing Murine Cardiomyocytes[END_REF]

(c n , b 1,n , b 2,n , b 3,n , b 4,n ) (L 2 (0,T ;H 1 (Ωcyt))∩L ∞ (0,T ;L 2 (Ωcyt))) 5 ≤ C, and (s n , b 5,n ) (L 2 (0,T ;H 1 (Ωsr))∩L ∞ (0,T ;L 2 (Ωsr))) 2 ≤ C,
for some constant C > 0 not depending on n. Finally, note that from (i) and (ii) we get easily for i = 1, . . . , 4

∂ t c n L 2 (0,T ;(H 1 (Ωcyt)) ) + ∂ t s n L 2 (0,T ;(H 1 (Ωsr)) ) + ∂ t b i,n L 2 (0,T ;(H 1 (Ωcyt)) ) + ∂ t b 5,n L 2 (0,T ;(H 1 (Ωsr)) ) ≤ C,
for some constant C > 0 not depending on n. Therefore, (iii) is a consequence of this, (ii) and the compactness argument (see for e.g. [START_REF] Simon | Compact sets in the space L p (0, T ; B)[END_REF] for more details). This concludes the proof of Lemma 2.5.

In summary, Lemma 2.5 and the theory of compact sets [START_REF] Simon | Compact sets in the space L p (0, T ; B)[END_REF], imply that there exists

(c, b 1 , b 2 , b 3 , b 4 , s, b 5 ) ∈ (L 2 (0, T ; H 1 (Ω cyt ))) 5 × (L 2 (0, T ; H 1 (Ω sr )
)) 2 such that, up to extracting subsequences if necessary,

(cn, b i,n , b 2,n , b 3,n , , b 4,n , b 5,n , sn) → (c, b 1 , b 2 , b 3 , b 4 , b 5 , s) in (L 2 (Ω cyt,T )) 5 × (L 2 (Ω sr,T )) 2
strongly and a.e., and in (L 2 (0, T ;

H 1 (Ωcyt))) 5 × (L 2 (0, T ; H 1 (Ωsr))) 2 weakly.
Consequently, the continuity of F on A holds. Moreover, using again Lemma 2.5 and the theory of compact sets [START_REF] Simon | Compact sets in the space L p (0, T ; B)[END_REF], we get the compactness of the map

F(A) → L 2 (Ω cyt,T ) × L 2 (Ω sr,T
). Then, thanks to Schauder fixed-point theorem, the operator F has a fixed point (c ε , s ε ). That is, there exists a solution to the approximate problem [START_REF] Louch | Control of Ca2+ Release by Action Potential Configuration in Normal and Failing Murine Cardiomyocytes[END_REF]. That is there exists a solution to

T 0 ∂ t c ε , ϕ c dt + Ω cyt,T D c (c ε )∇c ε • ∇ϕ c dx dt - T 0 Γryr g(s ε -c ε )ϕ c σdt = - Ω cyt,T 4 i=1 R i,ε (c ε , b i,ε )ϕ c dx dt, T 0 ∂ t s ε , ϕ s dt + Ω sr,T D s (s ε )∇s ε • ∇ϕ s dx dt + T 0 Γryr g(s ε -c ε )ϕ s σdt = - Ω sr,T R s,ε (s ε , b 5,ε )ϕ s dx dt, T 0 ∂ t b i,ε , ϕ i dt + D bi Ω cyt,T ∇b i,ε • ∇ϕ i dx dt = Ω cyt,T R i,ε (c ε , b i,ε )ϕ i dx dt, T 0 ∂ t b 5,ε , ϕ 5 dt + D b5 Ω sr,T ∇b 5,ε • ∇ϕ 5 dx dt = Ω sr,T R 5,ε (c ε , b 5,ε )ϕ 5 dx dt, (27)
for all ϕ c , ϕ i ∈ L 2 (0, T ; H 1 (Ω cyt )) and ϕ s , ϕ 5 ∈ L 2 (0, T ; H 1 (Ω sr )) for i = 1, . . . , 4.

Existence of weak solutions

We have shown in Section 2.1 that problem [START_REF] Pods | Electrodiffusion Models of Neurons and Extracellular Space Using the Poisson-Nernst-Planck Equations-Numerical Simulation of the Intra and Extracellular Potential for an Axon Model[END_REF]

admits a solution (c ε , b 1,ε , b 2,ε , b 3,ε , b 4,ε , s ε , b 5,ε
). The goal in this subsection is to send the regularization parameter ε to zero in sequences of such solutions to obtain weak solutions of the original system [START_REF] Lagnese | Dynamic Domain Decomposition in Approximate and Exact Boundary Control in Problems of Transmission for Wave Equations[END_REF]. Observe that, for each fixed ε > 0, we have shown the existence of a solution.

(c ε , b i,ε , b 2,ε , b 3,ε , b 4,ε , s ε , b 5,ε ) to (27) such that for i = 1, . . . , 4 0 ≤ c ε (x c , t), b i,ε (x c , t), s ε (x s , t), b 5,ε (x s , t) ≤ M,
for a.e. (x c , t) ∈ Ω cyt,T and (x s , t) ∈ Ω sr,T . Herein, the constant M > 0 is not depending on ε. Now, we substitute ϕ c = c ε , ϕ s = s ε , ϕ i = b i,ε in ( 27) for i = 1, . . . , 5, and we work exactly as in the proof of (ii) in Lemma 2.5, to obtain for i = 1, . . . , 4

sup 0≤t≤T Ω cyt |cε(x, t)| 2 dx + sup 0≤t≤T Ωsr |sε(x, t)| 2 dx + sup 0≤t≤T Ω cyt |bi,ε(x, t)| 2 dx + sup 0≤t≤T Ωsr |b5,ε(x, t)| 2 dx ≤ C, (28)
and

Ω cyt,T |∇cε(x, t)| 2 dx dt + Ω sr,T |∇sε(x, t)| 2 dx dt + Ω cyt,T |∇bi,ε(x, t)| 2 dx dt + Ω sr,T |∇b5,ε(x, t)| 2 dx dt ≤ C, (29)
for some constant C > 0 independent of ε. Repeating the steps of the proof of (iii) in Lemma 2.5, we derive the bound for i = 1, . . . , 4,

∂ t c ε L 2 (0,T ;(H 1 (Ωcyt)) ) + ∂ t s ε L 2 (0,T ;(H 1 (Ωsr)) ) + ∂ t b i,ε L 2 (0,T ;(H 1 (Ωcyt)) ) + ∂ t b 5,ε L 2 (0,T ;(H 1 (Ωsr)) ) ≤ C, (30)
for some constant C > 0. Then, combining ( 28), ( 29) and ( 30) with standard compactness results (cf. [START_REF] Simon | Compact sets in the space L p (0, T ; B)[END_REF]) we can extract subsequences, which we do not relabel, such that, as ε goes to 0, 2 , strongly and a.e., and in (L 2 (0, T ;

           (cε, b i,ε , sε, b 5,ε ) → (c, b i , s, b 5 ) in (L 2 (Ω cyt,T ) 2 ) × (L 2 (Ω sr,T ))
H 1 (Ωcyt))) 2 × (L 2 (0, T ; H 1 (Ωsr))) 2 weakly, (∂tcε, ∂tb i,ε ) → (∂tc, ∂tb i ) weakly in (L 2 (0, T ; (H 1 (Ωcyt)) )) 2 , (∂tsε, ∂tb 5,ε ) → (∂ts, ∂tb 5 ) weakly in (L 2 (0, T ; (H 1 (Ωsr)) )) 2 , (31)
for i = 1 . . .

which in turn implies that for i = 1, . . . , 5

R i,ε (c ε , b i,ε ) → R i (c, b i ), R s,ε (s ε , b 5,ε ) → R s (s, b 5 ), (33)
a.e. and strongly in L p for 1 ≤ p < ∞. As a consequence of ( 31), ( 32) and [START_REF] Vest | Defective cardiac ryanodine receptor regulation during atrial fibrillation[END_REF], it is clear that as ε → 0 in (27), we can identify the limit as the weak solution from Definition 2.2.

Remark 2. Note that, since the interaction terms are Lipschitz (recall that our solutions are uniformly bounded in L ∞ ) and the positivity of diffusion function, the uniqueness of weak solution can be obtained easily (for a similar proof see [START_REF] Anaya | On a vorticity-based formulation for reaction-diffusion-Brinkman systems[END_REF]).

Numerical approximation

In this section we construct our finite element scheme for the system [START_REF] Lagnese | Dynamic Domain Decomposition in Approximate and Exact Boundary Control in Problems of Transmission for Wave Equations[END_REF]. We establish several a priori (discrete energy) estimates for our scheme and compactness of discrete solutions which eventually will imply the desired convergence results to the unique weak solution.

The semi-discrete scheme

We consider a quasi-uniform family T i,h of triangulations of Ωi by tetrahedra K, for i = cyt, sr. The parameter h has the sense of an upper bound for the maximum diameter of the control volumes in T i,h for i = cyt, sr. For a given an integer k ≥ 0 and S ⊂ R 3 , we denote by P k (S) the space of polynomial functions defined in S of total degree up to k. Next, we define the finite element subspaces V i,h by

V i,h = m h ∈ H 1 (Ω i ) : m h | K ∈ P 1 (K) ∀K ∈ T i,h ,
where for any K ∈ T i,h (Ω i ) for i = cyt, sr. respectively. We consider the appropriate L 2 -projections of the initial condition (c 0 h , b 0 i,h ) = P V cyt,h (c 0 , b i,0 ) (for i = 1, . . . , 4) and (s 0 h , b 0 5,h ) = P V h (s 0 , b 5,0 ). To advance the numerical solution from t n-1 to t n = t n-1 + δt (δt = T /N ), we use the following implicit finite element scheme:

For n ∈ {1, . . . , N }, find (c n h , b n 1,h , b n 2,h , b n 3,h , b n 4,h ) ∈ (V cyt,h ) 5 and (s n h , b n 5,h ) ∈ (V sr,h ) 2 such that c n h -c n-1 h δt , ϕ c Ωcyt + Ωcyt D c (c n h )∇c n h • ∇ϕ c dx - Γryr g(s n h -c n h )ϕ c dσ = - 4 i=1 R i (c n h , b n i,h ), ϕ c Ωcyt , s n h -s n-1 h δt , ϕ s Ωsr + Ωsr D s (s n h )∇s n h • ∇ϕ s dx + Γryr g(s n h -c n h )ϕ s dσ = -R 5 (s n h , b n 5,h ), ϕ s Ωsr , b n i,h -b n-1 i,h δt , ϕ i Ωcyt + Ωcyt D bi ∇b n i,h • ∇ϕ i dx = R i (c n h , b n i,h ), ϕ i Ωcyt , b n 5,h -b n-1 5,h δt , ϕ 5 Ωsr + Ωsr D b5 ∇b n 5,h • ∇ϕ 5 dx = R 5 (s n h , b n 5,h ), ϕ 5 Ωsr , (34)
for all ϕ c , ϕ i ∈ V cyt,h and ϕ s , ϕ 5 ∈ V sr,h for i = 1, . . . , 4.

Numerical convergence analysis

Note that the problem of proving uniform boundedness of the discrete finite element solutions is still an open problem (recall that the positive and negative cuts test functions are not admissible in finite element context). So, comparing to the continuous case, we assume that the initial conditions (c ). We establish several a priori (discrete energy) estimates, which eventually will imply the desired convergence results.

Remark 3. In the proof of Lemmas 3.1 and 3.2 below, we will use the following Gagliardo-Nierenberg inequality:

For Ω ⊂ R 3 there exists a constant G such that for u ∈ H 1 (Ω)

u L q (Ω) ≤ G u α H 1 (Ω) u 1-α L 2 (Ω) for 1 q = 3 -2α 6 and α ∈ (0, 1). ( 35
) Lemma 3.1. Let (c n h , b n 1,h , b n 2,h , b n 3,h , b n 4,h , s n h , b n 5,h
) be part of the solution of [START_REF] Williams | Dynamics of Calcium Sparks and Calcium Leak in the Heart[END_REF]. Then there exists a constant C > 0 depending on c 0 L 2 (Ωcyt) , b i,0 L 2 (Ωcyt) , and s 0 L 2 (Ωsr) , b 5,0 L 2 (Ωsr) such that for i = 1, . . . , 4

(c h , b i,h ) L ∞ (0,T ;L 2 (Ωcyt,R 2 )) + (s h , b 5,h ) L ∞ (0,T ;L 2 (Ωcyt,R 2 )) ≤ C, ∇c h L 2 (Ω T ,cyt ) + 3 i=1 ∇b i,h L 2 (Ω T ,cyt ) + ∇s h L 2 (Ω T ,sr ) ≤ C. ( 36
)
Proof. We use [START_REF] Williams | Dynamics of Calcium Sparks and Calcium Leak in the Heart[END_REF] with

φ c = c n h , φ s = -s n h and φ i = b n i,h (for i = 1, . . . , 5
), and we sum over n = 1, . . . , k for all 1 < k ≤ N . The result is

1 2 (c k h , c k h ) Ω cyt + (s k h , s k h ) Ωsr + 4 i=1 (b k i,h , b k i,h ) Ω cyt + (b k 5,h , b k 5,h ) Ωsr + D min k n=1 Ω cyt ∇c n h 2 + k n=1 Ωsr ∇s n h 2 + Γryr g |s n h -c n h | 2 dσ ≤ 1 2 (c 0,h , c 0,h ) Ω cyt + (s 0,h , s 0,h ) Ωsr + 4 i=1 (b i,0,h , b i,0,h ) Ω cyt + (b 5,0,h , b 5,0,h ) Ωsr + C 1 k n=1 Ω cyt |c n h | 2 + 4 i=1 b n i,h 2 + k n=1 Ωsr |s n h | 2 + b n 5,h 2 + C 2 4 i=1 k n=1 Ω cyt (c n h ) 2 b n i,h + c n h (b n i,h ) 2 + k n=1 Ωsr (s n h ) 2 b n 5,h + s n h (b n 5,h ) 2 := 1 2 I 1 + C 1 I 2 + C 3 I 3 ,
for some constant C 1 , C 2 > 0. Herein, we have used the positivity of D c and D s , and the convexity inequality

a n (a n -a n-1 ) ≥ |a n | 2 2 - a n-1 2 2 .
Observe that since the initial conditions are bounded in L 2 , we get

I 1 ≤ C 3 ,
for some constant C 3 > 0. Regarding the term I 3 , we use Gagliardo-Nieremberg inequality (35) for q = 4, Poincaré inequality and Young inequality to deduce for i = 1, . . . ,

4 k n=1 Ω cyt (c n h ) 2 b n i,h ≤ k n=1 c n h L 2 (Ω cyt ) c n h L 4 (Ω cyt ) b n i,h L 4 (Ω cyt ) ≤ G k n=1 c n h 5/4 L 2 (Ω cyt ) c n h 3/4 H 1 (Ω cyt ) b n i,h 1/4 L 2 (Ω cyt ) b n i,h 3/4 H 1 (Ω cyt ) ≤ D min 2 k n=1 c n h H 1 (Ω cyt ) b n i,h H 1 (Ω cyt ) + C 4 k n=1 c n h 5 L 2 (Ω cyt ) b n i,h L 2 (Ω cyt ) ≤ D min 4 k n=1 c n h 2 H 1 (Ω cyt ) + b n i,h 2
H 1 (Ω cyt ) + C 5 k n=1 c n h 6 L 2 (Ω cyt ) + b n i,h 6
L 2 (Ω cyt ) ,
for some constants C 4 , C 5 > 0. Similarly, we get for for i = 1, . . . , 4

k n=1 Ωcyt c n h (b n i,h ) 2 ≤ D min 4 k n=1 c n h 2 H 1 (Ωcyt) + b n i,h 2 H 1 (Ωcyt) + C 6 k n=1 c n h 6 L 2 (Ωcyt) + b n i,h 6
L 2 (Ωcyt) , k n=1 Ωsr (s n h ) 2 b n 5,h ≤ D min 4 k n=1 s n h 2 H 1 (Ωsr) + b n 5,h 2 H 1 (Ωsr) + C 7 k n=1 s n h 6 L 2 (Ωsr) + b n 5,h 6
L 2 (Ωsr) ,
and

k n=1 Ωcyt s n h (b n 5,h ) 2 ≤ D min 4 k n=1 s n h 2 H 1 (Ωsr) + b n 5,h 2 H 1 (Ωsr) + C 8 k n=1 s n h 6 L 2 (Ωsr) + b n 5,h 6
L 2 (Ωsr) ,
for some constants C 6 , C 7 , C 8 > 0. Collecting the previous inequalities, we obtain

Ωcyt c k h 2 + 4 i=1 b k i,h 2 + Ωsr s k h 2 + b k 5,h 2 + D min 2 k n=1 Ωcyt ∇c n h 2 + k n=1 Ωsr ∇s n h 2 ≤ 2C 3 + 2C 1 k n=1 Ωcyt |c n h | 2 + 4 i=1 b n i,h 2 + Ωsr |s n h | 2 + b n 5,h 2 + C 9 k n=1 c n h 6 L 2 (Ωcyt) + 4 i=1 b n i,h 6 L 2 (Ωcyt) + k n=1 s n h 6 L 2 (Ωsr) + 4 i=1 b n 5,h 6 L 2 (Ωsr) , ≤ 2C 3 + C 10 k n=1 Ωcyt |c n h | 2 + 4 i=1 b n i,h 2 + Ωsr |s n h | 2 + b n 5,h 2 × 1 + Ωcyt |c n h | 2 + 4 i=1 b n i,h 2 + Ωsr |s n h | 2 + b n 5,h 2 2 , (37)
for some constants C 9 , C 10 > 0. Now we set

U := 1 + Ωcyt |c h | 2 + 4 i=1 |b i,h | 2 + Ωsr |s h | 2 + |b 5,h | 2 .
Observe that (37) implies

U k ≤ C 11 + C 10 k n=1 (U n ) 3 for every 1 < k ≤ N ,
and for some constant C 11 > 0. Therefore by the discrete Gronwall inequality [START_REF] Tao | A generalization of discrete Gronwall inequality and its application to weakly singular Volterra integral equation of the second kind[END_REF] to U , yields

Ωcyt |c h | 2 + 4 i=1 |b i,h | 2 + Ωsr |s h | 2 + |b 5,h | 2 ≤ U (t) ≤ Const( (c 0 , b i,0 , s 0 , b 5,0 ) ) L 2 (Ωcyt,R 2 )×L 2 (Ωsr,R 2 ) ,
for i = 1, . . . , 4. Therefore, we get the first estimate of (36). Using this and (37) we get the second part of estimate (36). This concludes the proof of Lemma 3.1.

Now, we consider

(c h , , b i,h ) and (s h , b 5,h ) the piecewise affine in t functions in (W 1,∞ ([0; T ]; V cyt,h )) 2 and in (W 1,∞ ([0; T ]; V sr,h )) 2
, respectively, and interpolating (c n h , b n i,h ) n=0,...,N ⊂ (V cyt,h ) 2 and (s n h , b n 5,h ) n=0,...,N ⊂ (V sr,h ) 2 at the points (n∆t) n=0,...,N (recall the interpolation scheme

n h = n n-1 h + t -t n-1 ∆t (n n h -n n-1 h ) for n = c, s, b i ) for i = 1, . . . , 4. Then we have                ∂ t c h = div (D c (c h )∇c h ) - 4 i=1 R i (c h , b i,h ), in Ω cyt,T , ∂ t b i,h = div (D bi ∇b i,h ) + R i (c h , b i,h ), i ∈ {1, . . . , 4}, in Ω cyt,T , ∂ t s h = div (D s (s h )∇s h ) -R 5 (s h , b 5,h ), in Ω sr,T , ∂ t b 5,h = div (D b5 ∇b 5,h ) + R 5 (s h , b 5,h ), in Ω sr,T . (38)
In the next lemma, we establish the relative compactness in L 2 of the sequences

(c h , b 1,h , b 2,h , b 3,h , b 4,h , s h , b 5,h ),
which is achieved by constructing space and time translates and using the a priori estimates given above.

Lemma 3.2. There exist a constant C > 0 depending on initial conditions, Ω sr , Ω cyt and T such that for i = 1, . . . , 4

(c h , b i,h ) (L 2 (0,T ;H 1 (Ωcyt))) 2 + (s h , b 5,h ) (L 2 (0,T ;H 1 (Ωsr))) 2 ≤ C, (39)
and

T 0 ∂ t c h , ϕ c cyt dt + 4 i=1 T 0 ∂ t b i,h , ϕ i cyt dt + T 0 ∂ t s h , ϕ s sr dt + T 0 ∂ t b 5,h , ϕ 5 sr dt ≤ C ϕ c L 2 (0,T ;H 1 (Ωcyt)) + 4 i=1 ϕ i L 2 (0,T ;H 1 (Ωcyt)) + ϕ s L 2 (0,T ;H 1 (Ωsr)) + ϕ 5 L 2 (0,T ;H 1 (Ωsr)) , (40)
for any ϕ c , ϕ i ∈ L 2 (0, T ; H 1 (Ω cyt )) and ϕ s , ϕ 5 ∈ L 2 (0, T ; H 1 (Ω sr )).

Proof. First, we begin by proving (39) but this is the consequence of the definition of (c h , s h , b i,h ) (for i = 1, . . . , 5) and Lemma 3.1.

The next part of the prove is devoted to (40). Let ϕ c ∈ L 2 (0, T ; H 1 (Ω)). Using (38), Lemma 3.1 Gagliardo-Nieremberg inequality, Poincaré inequality and Young inequality, we obtain

T 0 ∂ t c h , ϕ c cyt dt = T 0 (D c (c h )∇c h , ϕ c ) cyt dt - 4 i=1 T 0 (R i (c h , b i,h ), ϕ c ) cyt dt ≤ D max ∇c h L 2 (Ω cyt,T ) ∇ϕ c L 2 (Ω cyt,T ) + C 12 c h L 2 (Ω cyt,T ) + 4 i=1 b i,h L 2 (Ω cyt,T ) ϕ c L 2 (Ω cyt,T ) + C 13 4 i=1 T 0 (|c h , b i,h | , |ϕ c |) cyt dt ≤ C 14 ϕ c L 2 (0,T ;H 1 (Ωcyt)) + C 13 4 i=1 I i , (41)
for some constants C 12 , C 13 , C 14 > 0. Regarding the integral I i , we use Gagliardo-Nieremberg inequality (35) (for q = 3) and Young inequality to obtain for i = 1, . . . , 4

I i ≤ T 0 c h L 2 (Ωcyt) b i,h L 3 (Ωcyt) ϕ c L 6 (Ωcyt) dt ≤ G c h L ∞ (0,T ;L 2 (Ωcyt)) b i,h L ∞ (0,T ;L 2 (Ωcyt)) T 0 b i,h H 1 (Ωcyt) ϕ c L 6 (Ωcyt) dt ≤ G c h L ∞ (0,T ;L 2 (Ωcyt)) b i,h L ∞ (0,T ;L 2 (Ωcyt)) × b i,h L 1 (0,T ;H 1 (Ωcyt) ϕ c L 2 (0,T ;L 6 (Ωcyt)) dt. (42)
Now we use Lemma 3.1 and the continuous embedding of L 2 (0, T ; H 1 (Ω cyt )) into L 2 (0, T ; L 6 (Ω cyt )) to get from (42)

I i ≤ C ϕ c L 2 (0,T ;H 1 (Ωcyt)) ,
for some constant C > 0, for i = 1, . . . , 4. From this and (41) we deduce the bound

∂ t c h L 2 (0,T ;(H 1 (Ω)) ) ≤ C. (43)
Reasoning similarly, we get (40) for (b i,h , s h ) for i = 1, . . . , 5, therefore finishing our proof.

Observe that from the definition of (c h , s h , b i,h ) together with [START_REF] Williams | Dynamics of Calcium Sparks and Calcium Leak in the Heart[END_REF] eventually implies that

c h -c h 2 L 2 (Ω cyt,T ) + 4 i=1 b i,h -b i,h 2 L 2 (Ω cyt,T ) ≤ N n=1 ∆t c n h -c n-1 h 2 L 2 (Ωcyt) + 4 i=1 b n i,h -b n-1 i,h 2
L 2 (Ωcyt) ≤ C(∆t) → 0 as ∆t → 0,
and

s h -s h 2 L 2 (Ω sr,T ) + b 5,h -b 5,h 2 L 2 (Ω sr,T ) ≤ N n=1 ∆t s n h -s n-1 h 2 L 2 (Ωsr) + b n 5,h -b n-1 5,h 2 L 2 (Ωsr) ≤ C(∆t) → 0 as ∆t → 0.
Therefore, from Lemma 3.2 and standard compactness results (see for e.g. [START_REF] Simon | Compact sets in the space L p (0, T ; B)[END_REF]), there exists a subsequence of

(c h , s h , b 1,h , b 2,h , b 3,h , b 4,h , b 5,h ) not relabeled, such that, as h → 0,          (c h , b 1,h , b 2,h , b 3,h , b 4,h ) → (c,

Numerical simulation

In this section, we present some numerical simulations in several experiments with different parameters. In the first experiment, we investigate the diffusion function to catch the accurate calcium wave front characteristics using imaging observations as reference [START_REF] Chen | Calcium waves initiating from the anomalous subdiffusive calcium sparks[END_REF]. In the next experiment, we analyze the influence of the RyR conductance on calcium front wave propagation in cytosolic domain, followed by an investigation on pathological consequences on the cellular scale. Finally, we extend our simulation results in three dimensional framework and we discuss the phenomenological behavior of calcium spark through RyR.

To approximate the solution, we consider the finite element framework defined in Section 3 along with Euler implicit time discretization with time step δt = 10 -5 s. The mesh parameters are determined using mesh adaptation to deal with high gradient of the solutions that came from the non uniform distribution of b 4 and b 5 (See Figure 3b below). For calcium initial concentration, we take c 0 = 0.14 in Ω cyt and s 0 = 1.3 × 10 3 in Ω sr .

For buffers initial concentration in Ω cyt we put b 1,0 = 318 × 10 -3 , b 2,0 = 471 × 10 -3 , b 3,0 = 2.82, b 4,0 = 13.2 and the initial concentration of CSQN in Ω sr is b 5,0 = 24.6 × 10 3 . The values of previous parameters are provided in µM . We present the diffusion coefficient and interaction parameters (see [START_REF] Hake | Modelling cardiac calcium sparks in a three-dimensional reconstruction of a calcium release unit[END_REF] and the references therein) in Table 1. Almost all parameters have been discovered in cardiac cell medium. 1 shows experimental data related to cardiac cell. The units are presented as follows : The diffusion coefficient (D l ) unit is 10 3 nm 2 ms -1 and the total mass unit (B i ) is µM . The k on i unit is 10 -3 ms -1 µm -1 and k of f i unit is ms -1 . The diffusion of buffers b 4 and b 5 are taken very small.

2D numerical simulation of calcium model

In cell medium, we mention that fixed buffers are not uniformly distributed (see [START_REF] Hake | Modelling cardiac calcium sparks in a three-dimensional reconstruction of a calcium release unit[END_REF] for more details) which would induce a higher gradient of the solution in some regions. Hence, we consider a random distribution of buffers and we adapt the mesh generation (see Figure 3) to achieve a finer approximation in the critical regions. The mesh adaptation is performed using a maximal edge length h max = 0.3 on both computational domain with maximal number of vertexes 120 × 10 3 (resp. 200 × 10 3 ) for cytosolic domain (resp. sarcoplasmic domain). Figure 4 shows the two dimensional Ca 2+ -distribution for a time instant of t = 1ms after the calcium spark. The concentration dependent diffusion coefficient is given by :

D s (v) = 0.1(1 + v) and D c (u) = D c (0.5 + u m )
with different values of m (A similar formulation is given in [START_REF] Crank | The Mathematics of Diffusion[END_REF]). By changing m values from 0 to 5 (see Figure 5b), we observe that Ca 2+ wavefront displays an ellipsoidal shape with a fast decaying near to the front wave. In imaging results from [START_REF] Chen | Calcium waves initiating from the anomalous subdiffusive calcium sparks[END_REF], we observe the same phenomenon due to the heavy fixed buffering in the cell medium that decreases the calcium diffusion in the low concentration areas. We mention also the explosive increasing of release flux at the beginning of calcium spark which has been confirmed experimentally in [START_REF] Williams | Dynamics of Calcium Sparks and Calcium Leak in the Heart[END_REF]. Hence, The non constant diffusion coefficient induces the same phenomenon. The larger m is taken more the explosive increasing is important. To understand the wavefront differences between linear (D c (u) = 1.5D c and D s (v) = D s ) and nonlinear case (m = 3), Figure 6 shows the evolution Ca 2+ -distribution in both cases. In the second row, we observe a large propagating front wave moving away from the Ca 2+ highest concentration around RyR (similar observation shown in [START_REF] Hake | Modelling cardiac calcium sparks in a three-dimensional reconstruction of a calcium release unit[END_REF]) because of the linear diffusion with the simple Fickian law. On the upside, the propagating front wave is much smaller and clear in this case. Moreover, this affects the release flux through RyR, see Figure 5a. Figure 7, shows the influence of RyR conductance on calcium wave front velocity. Figure 7c shows a faster wave front than in Figure 7a case. Pathological behaviors of RyR influences directly the calcium wave velocity and consequently influences the heart beat rate. In other cases some perturbations on heavy buffering protein (CSQN) inside the SR influence the release flux through RyR. In this case, the universal RyR conductance should be corrected to regularize the calcium wave by regularizing the release flux through ionic channels.

3D numerical simulation of calcium model

To fulfill the requirement of the finite element method in three dimensional space we consider a tetrahedral mesh of both domain. Figure 8 shows a section of the external domain mesh (Cytosolic domain) and a meshing of the sarcoplasmic domain. The right figure present the section of Ω sr where the visualization of Ca 2+ -concentration occur.

Figure 8. In the left hand-side we show complete mesh grid of Ω sr and a section of mesh grid of Ω cyt . In the right hand-side we show the volume section in which 3D simulation will be visualized.

In Figure 9, we present a numerical visualization of calcium solution in both domains. In the sarcoplasmic domain, we show the calcium distribution on the internal sarcoplasmic membrane. We choose to present the calcium dynamics in two different 2D sections with different distance from the ionic channel. To clarify the calcium wave sharpness, we have depict different isosurfaces. Each one describes a specific value of calcium concentration. The value of isosurfaces equidistantly varying between 0.14 and 0.44 for linear and nonlinear diffusion cases. In Figure 10, we plot the discs shown in Figure 9. The white curves (in Figure 10) present the line curves attached to values of calcium concentrations. To understand the calcium wavefront differences between the anomalous diffusion model and a linear model, we have computed the calcium distribution in both cases in different time frames. Here, we clearly notice the main difference between the linear diffusion and nonlinear diffusion. In Figure 9, the isosurfaces present a splaying dynamics due to the linear diffusion in the medium comparing with the anomalous diffusion. The spatial dynamics in Figure 10 clarifies the anomalous diffusion between the linear model (see for e.g. [START_REF] Chai | Towards simulation of subcellular calcium dynamics at nanometre resolution[END_REF]) and the proposed model. In the linear model we observe that calcium concentration is much higher near to the ionic channels than further (to the ionic channels) which has been shown in [START_REF] Chai | Towards simulation of subcellular calcium dynamics at nanometre resolution[END_REF]. But this doesn't cover the heterogeneity of the cell medium since the constant diffusion based on the assumption that the calcium flux is constantly proportional to the calcium gradient in the whole domain. In the right two columns of Figure 10, we observe that the calcium front-wave is sharper than the linear case. The choice of an increasing function of calcium concentration may explain the slow penetration of calcium through the cellular heterogeneous media (we refer to [START_REF] Gal | Experimental evidence of strong anomalous diffusion in living cells[END_REF][START_REF] Severs | Remodelling of gap junctions and connexin expression in diseased myocardium[END_REF] for more details). The complex dense cell structure prevents normal diffusion in Ca 2+ low concentration, and the fast diffusion occurs in high concentration region due to high kinetics of Ca 2+ particles (for more information see [START_REF] Crank | The Mathematics of Diffusion[END_REF] and the references therein), which should form a sharp wavefront pattern. This is approved by experimental micro-imaging technics on calcium diffusion in the cardiac myocyte.

Finally, to test the numerical scheme convergence, we consider a simple geometry on Ω sr = [-2, 0] × [0, 1] and Ω cyt = [0, 4] × [0, 1]. In space triangularization, we consider a uniform mesh grid of length h = 1 N with N varying from 5 to 250, the time step δt = 0.05 and the final time value T = 2. Table 2 shows errors relative to the exact solutions using P 1 and P 2 finite element. To simplify, we do not consider the interaction with buffers. We have used as an exact solution (for more information, see [START_REF] Jaeger | An Investigation of Necessary Grid Resolution for Numerical Simulations of Calcium Dynamics in Cardiac Cells[END_REF]) the following function : After running the test using different parameters. We present in Table 2

Conclusion & Perspective

In this work, we have proposed a nonlinear reaction diffusion model describing the evolution of calcium in cardiac cell. Based on compactness results, we have studied the well posedness of the problem. Then, we have used a finite element scheme along with implicit Euler discretization to approximate the biochemical substance concentrations. We attached a convergence analysis of our scheme using Kolmogorov compactness results. Moreover, we have produced numerical simulations to clarify the physiological link between the given parameters and calcium behavior under linear and anomalous diffusion. On the experimental scale, B. S. Donahue and R. F. Abercrombie [START_REF] Donahue | Free diffusion coefficient of ionic calcium in cytoplasm[END_REF] deduced that the calcium diffusion in cytosol depends on calcium concentration. Advanced micro-imaging technology in cardiology provides sophisticated results on calcium diffusion profile. Figure 6 proves encouraging similarities between Ca 2+ -distribution and imaging results shown in [START_REF] Chen | Calcium waves initiating from the anomalous subdiffusive calcium sparks[END_REF]. A precise prediction of spacial distribution (including nonlinear diffusion properties) in cytosol is widely important, it may offer some explanation of spontaneous Ca 2+ -sparks since passive ionic receptor are very sensitive to local calcium concentration. The linear models [START_REF] Chai | Towards simulation of subcellular calcium dynamics at nanometre resolution[END_REF][START_REF] Hake | Modelling cardiac calcium sparks in a three-dimensional reconstruction of a calcium release unit[END_REF][START_REF] Langer | Calcium Concentration and Movement in the Diadic Cleft Space of the Cardiac Ventricular Cell[END_REF][START_REF] Louch | Control of Ca2+ Release by Action Potential Configuration in Normal and Failing Murine Cardiomyocytes[END_REF][START_REF] Peskoff | Calcium Concentration and Movement in the Ventricular Cardiac Cell during an Excitation-Contraction Cycle[END_REF][START_REF] Peskoff | sarcolemmal calcium binding sites in heart: II. Mathematical model for diffusion of calcium released from the sarcoplasmic reticulum into the diadic region[END_REF] has shown some paradoxical results with respect to experiment, our model offers a general framework that may give a step ahead to solve incompatibilities between classical models and experiment. In further work, we intend to simulate the therapeutic strategy by targeting the ionic channel's characteristics and establishing a nonlinear optimal control problem on calcium dynamics.

Figure 1 .

 1 Figure 1. Description of calcium cycle in Cardiac cell

  c). Due to Ficks law, we have D s (s)∇s • η s = g(s -c), on Γ ryr,T := Γ ryr × (0, T ), and D c (c)∇c • η c = -g(s -c) on Γ ryr,T := Γ ryr × (0, T ),

Figure 2 .

 2 Figure2. A Graphical representations of the boundary of the calcium model. We note the sarcoplasmic membrane Γ sr = Γ i ∪ Γ ryr and cytosolic membrane Γ = Γ e ∪ Γ i ∪ Γ ryr , while we consider Γ cyt iso = Γ e ∪ Γ i the impermeable part of the cytosolic membrane, Γ sr iso = Γ i is the impermeable part of the sarcoplasmic membrane.

, 4 . 1 )

 41 Here, •, • denotes the duality pairing between H 1 and (H Definition 2.2 (L ∞ bound solution). A weak L ∞ bound solution of (16) is a nonnegative and a L ∞ function with seven components (c, b 1 , b 2 , b 3 , b 4 , s, b 5 ) satisfying Definition 2.1.

Theorem 2 . 3 .

 23 Assume that condition (17) holds. If c 0 , b 1,0 , b 2,0 , b 3,0 , b 4,0 ∈ L ∞ (Ω cyt ) and s 0 , b 5,0 ∈ L ∞ (Ω sr ), then the system (16) possesses a unique weak solution in sense of Definition 2.2.

Lemma 2 . 4 .

 24 Let c 0 , b 1,0 , b 2,0 , b 3,0 , b 4,0 ∈ L ∞ (Ω cyt ) and s 0 , b 5,0 ∈ L ∞ (Ω sr ), then for any (ĉ, ŝ) ∈ L 2 (Ω cyt,T ) × L 2 (Ω sr,T ), the system (19) has a unique solution. Moreover, there exists a constant C > 0 depending on the initial conditions, such that

  by c n , b i,n (for i = 1, . . . , 4), s n and b 5,n respectively. Herein, we replace the solution (c, b 1 , b 2 , b 3 , b 4 , s, b 5 ) by (c n , b 1,n , b 2,n , b 3,n , b 4,n , s n , b 5,n ). Next, we sum and integrate the resulting equation on (0, T ), on Ω cyt for (c n , b 1,n , b 2,n , b 3,n , b 4,n ) and on Ω sr for (s n , b 5,n ). Hence the Gronwall lemma and Poincare inequality yields

, 4 .

 4 This result, together with the weak-convergences inL ∞ of (c ε , b 1,ε , b 2,ε , b 3,ε , s ε ) to (c, b 1 , b 2 , b 3 , s), yields (cε, b 1,ε , b 2,ε , b 3,ε , b 4,ε ) → (c, b 1 , b 2 , b 3 , b 4 ), (sε, b 5,ε ) → (s, b 5) strongly in (L p (Ω cyt,T ))5 , and in (L p (Ω sr,T )) 2 , respectively, for 1 ≤ p < ∞,

  0 , b 1,0 , b 2,0 , b 3,0 , , b 4,0 ) and (s 0 , b 5,0 ) are only bounded in L 2 . The convergence of solutions generated by the scheme (34) is next established following two main lemmas. The first one states L ∞ (L 2 ) ∩ L 2 (H 1 ) boundedness of the discrete concentrations and buffers (c h , b 1,h , b 2,h , b 3,h , b 4,h , s h , b 5,h

  b 1 , b 2 , b 3 , b 4 ) strongly and a.e. in (L 2 (Ω cyt,T ) 5 , (s, b 5,h ) → (s, b 5 ) strongly and a.e. in L 2(Ω sr,T ) 2 , (c h , b 1,h , b 2,h , b 3,h , b 4,h ) → (c, b 1 , b 2 , b 3 , b 4 ) weakly in L 2 (0, T ; H 1 (Ω cyt )) 5 , (s, b 5,h ) → (s, b 5 ) weakly in L 2 (0, T ; H 1 (Ω sr )) 2 .(44) Finally, from the convergence in (44) we can identify the limit (c, b 1 , b 2 , b 3 , b 4 , s, b 5 ) as the weak solution of (16), and the convergence result is summarized as follows. Theorem 3.3. Assume that condition (17) holds. If c 0 , b 1,0 , b 2,0 , b 3,0 , b 4,0 ∈ L 2 (Ω cyt ) and s 0 , b 5,0 ∈ L 2 (Ω sr ), then the finite element solution (c h , b 1,h , b 2,h , b 3,h , b 4,h , s h , b 5,h ), generated by (34), converges along a subsequence to (c, b 1 , b 2 , b 3 , b 4 , s, b 5 ) as h, δt → 0, where (c, b 1 , b 2 , b 3 , b 4 , s, b 5 ) is a weak solution of (16) in the sense of Definition 2.1.

Figure 3 .

 3 Figure 3. The left hand-side figure represents the distribution of b 4 and b 5 . The right hand-side figure shows the Meshing used during the simulation

Figure 4 .

 4 Figure 4. The Graphical representation of the approximated solution in Cytosolic domain with non constant diffusion coefficient with m = 0, ..., 5 at t = 50 × 10 -5 s.

Figure 5 .

 5 Figure 5. Figure 5a shows temporal dynamics of release flux (i.e. Γryr g(s -c)) and Figure 5b represents a one dimensional calcium spatial distribution for various values of m = 1, ..., 5

Figure 6 .

 6 Figure 6. The figures on the down side represent the evolution of calcium distribution in the linear case. The figures on the up side devoted to the nonlinear case with m = 3

Figure 7 .

 7 Figure 7. The influence of RyR conductance on calcium wave front velocity. The values of RyR conductance g in Figures 7a, 7b and 7c are respectively 50, 150 and 300

5 Figure 9 .

 59 Figure 9. 3D presentation of calcium distribution in both domain.

Figure 10 .

 10 Figure 10. The first (resp. second) left column present the linear Ca 2+ distribution near (resp. far from) of disc to the ionic channel. The third (resp. forth) column presents the nonlinear distribution of Ca 2+ near (resp. far from) of the ionic channel.

8 = k 1 π 4 , b 8 = k 1 π 2 , 2 = k s a 2 3 , b 2 = a 2 , b 3 =D s a 2 3 D c , b 4 = a 4

 848223223344 s e (t, x, y) = b 1 exp(-b 2 t)sin(b 3 (x + 2)) + b 4 x + b 5 + b 6 cosh(b 7 x)cos(b 8 y)c e (t, x, y) = a 1 exp(-a 2 t) sin(a 3 (x -4)) + a 4 x + a 5 + a 6 cosh(a 7 x) cos(a 8 y)The functions c e (t, x, y) and s e (t, x, y) satisfies the transmission boundary problem for a specific constant. Let A = 100, C = -10, k 1 = 2, k 2 = 2, a 1 = -310, a 3 = -0.0072, a 4 = -1.97. The other constants are defined as follows :a a7 = a 8 , b 7 = b 8 , a 6 = A cosh(4a 7 ) , a 6 = C cosh(2b 7 ) , a D s D c , b 1 = D s a 1 a 3 cos(4a 3 ) D c b 3 cos(2b3 ) and a 5 = s 0 + 4a 4 , b 5 = c 0 -2b 4 For the conductance value: g(t, y) = D s a 1 a 3 exp(-a 2 t) cos(4a 3 ) + a 4 b 1 exp(-b 2 t) sin(-2b 3 ) + b 5 + b 6 cos(b 3 y)

Table 1 .

 1 Total concentration, reaction and diffusion rates in a cardiac cell.

	Diffusion coefficient D c	D s	D b1	D b2	D b3
	Value	220 73.3	140	25	42
	Total concentration B 1	B 2	B 3	B 4	B 5
	Value	140 25	42	42	42
	Off rates	k of f 1	k of f 2	k of f 3	k of f 4	k of f 5
	Value	45	0.238 0.110 0.0196 0.065
	On rates	k on 1	k on 2	k on 3	k on 4	k on 5
	Value	255 34	110	32.7	102

Table 2 .

 2 different values of the errore 2 L 2 = Ω cyt,T |c e -c h | 2 dx dt + Ω sr,T |s e -s h | 2 dx dt |c e -c h | 2 dx dt+ Ω cyt,T |∇c e -∇c h | 2 dx dt+ Ω sr,T |s e -s h | 2 dx dt+ Ω sr,T |∇s e -∇s h | 2 dx dt. This table represents an estimate of the error between the exact solution and the approximated solution

	and
	e 2 H 1 =
	Ω cyt,T
	Hence, we have the following table :
		P 1		P 2
	N	e L 2	e H 1	e L 2	e H 1
	20	0.163607	1.20455	0.00106206	0.0910758
	40 0.0413642	0.307202	6.92604E-5 0.0121671
	60 0.0184233	0.137055	1.3835E-5 0.00368335
	100 0.00663965 0.049437	1.80705E-6 0.000809208
	160 0.00259459 0.0193244 2.76497E-7 0.000199436
	250 0.0010629 0.00791734 5.09595E-8 5.25769E-5

research project for supporting this work.


domain: info math

text:
Miguel Couceiro

AND Lucien Haddad

Karsten Sch Ölzel

ON THE LOWER PART OF THE LATTICE OF PARTIAL CLONES

Let k be a k-element set. We show that the lattice of all strong partial clones on k has no minimal elements. Moreover, we show that if C is a strong partial clone, then the family of all partial subclones of C is of continuum cardinality. Finally we show that every non-trivial strong partial clone contains a family of continuum cardinality of strong partial subclones.

To the memory of Professor Ivan Stojmenovic

Preliminaries

In this section we recall basic concepts and terminology needed throughout the paper. For further background, see [START_REF] Börner | Maximal Partial Clones with no finite basis[END_REF][START_REF] Couceiro | A Survey on Intersections of Maximal Partial Clones of Boolean Partial Functions[END_REF][START_REF] Haddad | Countable Intervals of Partial Clones[END_REF][START_REF] Lau | Function algebras on finite sets. A basic course on many-valued logic and clone theory[END_REF].

Let A be a finite set with |A| ≥ 2. If |A| = k ≥ 2, then we assume that A = k := {0, . . . , k -1}. For a positive integer n, an n-ary partial function on k is a map f : dom(f ) → k where dom(f ) is a subset of k n called the domain of f . We denote by Par (n) For positive integers n, m ≥ 1, and f ∈ Par (n) (k) and g 1 , . . . , g n ∈ Par (m) (k), the composition of f with g 1 , . . . , g n is the partial function f [g 1 , . . . , g n ] ∈ Par (m) (k) whose domain dom(f [g 1 , . . . , g n ]) is the set of all a ∈ m i=1 dom(g i ) such that (g 1 ( a), . . . , g m ( a)) ∈ dom(f ), and whose values are given by f [g 1 , . . . , g n ]( a) := f (g 1 ( a), . . . , g n ( a)) for all a ∈ dom(f [g 1 , . . . , g n ]). For every positive integer n and each 1 ≤ i ≤ n, we denote by e n i the n-ary i-th projection function (or, simply, projection) defined by e n i (a 1 , . . . , a n ) = a i for all (a 1 , . . . , a n ) ∈ k n . Furthermore, we use the notation J k := {e n i : 1 ≤ i ≤ n, n ∈ N} for the set of all projections. Definition 1. A partial clone on k is a subset of Par(k) closed under composition and containing the set of all projection functions J k .

Remark 2. Partial clones can be equivalently defined in the formalism of Mal'cev; see, e.g., Chapter 20 of [START_REF] Lau | Function algebras on finite sets. A basic course on many-valued logic and clone theory[END_REF].

The set of partial clones on k forms a lattice L P k under inclusion, in which the infimum is the set-theoretical intersection. Thus the intersection of arbitrarily many partial clones on k is itself a partial clone on k. Clearly, the set J k of all projections on k is the smallest element in the lattice L P k of partial clones.

For f, g ∈ Par (n) (k), the partial function g is said to be a subfunction of f , fact that we denote by g ≤ f , if

g = f | dom(g) ,
i.e., if dom(g) ⊆ dom(f ) and g( a) = f ( a) for all a ∈ dom(g). In what follows, by an n-ary partial projection we mean a subfunction of a projection function whose domain is a proper subset of k n .

Example 4. It is easy to see that the set {f ∈ Par(k) : (0, . . . , 0) ∈ dom(f ) =⇒ f (0, . . . , 0) = 0} is a strong partial clone on k.

Strong partial clones have been widely studied in the literature, see, e.g., [START_REF] Börner | Maximal Partial Clones with no finite basis[END_REF][START_REF] Couceiro | A Survey on Intersections of Maximal Partial Clones of Boolean Partial Functions[END_REF][START_REF] Haddad | Countable Intervals of Partial Clones[END_REF]] and [START_REF] Lau | Function algebras on finite sets. A basic course on many-valued logic and clone theory[END_REF] for a long list of references. The recent paper [START_REF] Schölzel | Dichotomy on intervals of strong partial Boolean clones[END_REF] is dedicated to the study of strong partial clones of Boolean functions. It is well known that a partial clone C is strong if and only if Str(J k ) ⊆ C (see, e.g., [START_REF] Börner | Maximal Partial Clones with no finite basis[END_REF], Lemma 2.11). Moreover, it is easy to see that the intersection of an arbitrary family of strong partial clones on k is also a strong partial clone. In fact, the set of all strong partial clones, ordered by inclusion, also forms a lattice L Str(P k ) whose smallest element is Str(J k ).

For F ⊆ Par(k), let pclone(F ) denote the strong partial clone generated by F , i.e., the intersection of all strong partial clones on k containing F . As it is well-known, strong partial clones coincide exactly with the sets of polymorphisms of relations on k, that we now briefly describe.

For m ≥ 1, let ρ be an m-ary relation on k (i.e., a subset of k m ).

Denote by M(ρ, dom(f )) the set of all m × n matrices M on k whose columns M * j ∈ ρ, for j = 1, . . . , n and whose rows M i * ∈ dom(f )

for i = 1, . . . , m. A partial function f ∈ Par (n) (k) is said to preserve ρ ⊆ k m if (f (M 1 * ), . . . , f (M m * )) ∈ ρ for all M ∈ M(ρ, dom(f )). Set pPol ρ := {f ∈ Par(k) : f preserves ρ}.
It is well known that pPol ρ is a strong partial clone on k called the partial clone determined by ρ. Note that f ∈ pPolρ if there is no matrix M with columns in ρ and rows in dom(f ).

Main achievements:

The main contributions of this paper are 3-fold.

Minimal Strong Partial Clones

A minimal partial clone on k is an atom of the lattice L P k . In other words, a partial clone C is minimal if there is no partial clone C 0 such that J k = C 0 ⊂ C ⊂ Par(k). Minimal partial clones consisting of partial functions have been completely described in [START_REF] Börner | Minimal Partial Clones[END_REF], none of which is strong. The natural question is then to ask whether strong minimal partial clones exist, i.e., whether atoms of L Str(P k ) exist.

In this section we show that such partial clones do not exist. This result is not surprising, but to our knowledge it has never been published.

The idea behind the proof we present here comes from a discussion between the second author and Ferdinand Börner just after the paper [START_REF] Börner | Minimal Partial Clones[END_REF] has appeared. We will show that, given a partial function g that is not a partial projection function, we can construct a partial function h ∈ pclone(g), with h ∈ Str(J k ) and such that g ∈ pclone(h).

Let g ∈ Str(J k ) be an n-ary partial function and C := pclone(g). We distinguish two cases:

Case 1. For every (a 1 , . . . , a n ) ∈ dom(g), g(a 1 , . . . , a n ) ∈ {a 1 , . . . , a n }.

Let f be a restriction of g so that any further restriction of the domain of f results in a partial projection. Let m := |dom(f )| and let α 1 , . . . , α m ∈ dom(f ) ⊆ k n with α i := (a i 1 , . . . , a i n ) and such that

f      a 1 1 . . . a 1 n a 2 1 . . . a 2 n . . . . . . . . . a m 1 . . . a m n      =     b 1 b 2 . . . b m     . ( 1
)
Since f is not a partial projection, the right tuple (b 1 , . . . , b m ) T is not a column of the matrix on the left. Let M be the m × n matrix in [START_REF] Börner | Minimal Partial Clones[END_REF] and let ρ M be the m-ary relation on k consisting of all columns of M . Clearly f ∈ pPolρ M , and thus g ∈ pPolρ M .

Let f be the (n + 2)-ary partial function given by

f (x 1 , . . . , x n+2 ) := f (x 1 , . . . , x n ) for all (x 1 , . . . , x n+2 ) ∈ dom(f ) × k 2 . Note that f ∈ pclone(f ) ⊆ pclone(g).
Now take a c ∈ k, c = b 1 , and consider the subfunction h of f whose domain is the set of all rows of the following matrix:

       a 1 1 . . . a 1 n c b 1 a 1 1 . . . a 1 n b 1 c a 2 1 . . .
       =       b 1 b 1 b 2 . . . b m       . Note that |dom(h)| = m + 1.
We show that any restriction of h to a set with at most m tuples is a partial projection. Let h be a restriction of h with |dom(h )| ≤ m.

If (a 1 1 , . . . , a 1 n , c, b 1 ) ∈ dom(h ), then h ≤ e n+2 n+1 , i.e., h is a partial pro- jection on its (n + 1)-st coordinate. Similarly, if (a 1 1 , . . . , a 1 n , b 1 , c) ∈ dom(h ), then h ≤ e n+2
n+2 . On the other hand if (a i 1 , . . . , a i n , b i , b i ) ∈ dom(h ) for some i = 2, . . . , m, then h is a partial projection function on the same coordinate as the subfunction of f whose domain is contained in the set of rows of M except the i-th row. Now we show that h ∈ pPolρ M . Let N be an m × (n + 2) matrix with all columns in ρ M and all rows in dom(h). Then as N has m rows, the restriction of h on the rows of N is a partial projection function and so it preserves any relation, including the relation

ρ M . Thus h ∈ pPolρ M , i.e., pclone(h) ⊆ pPolρ M . Since g ∈ pPolρ M we conclude that g ∈ pclone(h) and thus pclone(h) C. Case 2. There is an n-tuple (a 1 , . . . , a n ) ∈ dom(g) such that g(a 1 , . . . , a n ) ∈ {a 1 , . . . , a n }.

Clones of Partial Projections

In this section we consider partial clones consisting only of projections and partial projection functions. For F ⊆ Par(k), let F denote the partial clone generated by F , i.e., the intersection of all partial clones on k containing F . Notice that F is not necessarily a strong partial clone.

Let R n denote the set of all n-tuples on k of weight 1, i.e.,

R n := {(x 1 , . . . , x n ) ∈ k n : n i=1 x i = 1}.
Furthermore, let p n be the n-ary partial projection given by

p n (x 1 , . . . , x n ) := x 1 if (x 1 , . . . , x n ) ∈ R n undefined otherwise.
Let P := {p n : n ≥ 2}. We will show that for t ≥ 2, p t ∈ P \ {p t } .

Recall Lemma 2.10 in [START_REF] Börner | Maximal Partial Clones with no finite basis[END_REF].

Lemma 6. Let F ⊆ Par(k) and set D 0 := F ∪ J k . For ≥ 0, define D +1 := D ∪ {f [g 1 , . . . , g n ] | f ∈ D (n) 0 and g 1 , . . . , g n ∈ D (m)
, for some n, m ≥ 1}.

Then F = ≥0 D .
Fix t ≥ 2 and F := P \ {p t }. Observe that, since F ⊆ Str(J k ), F ⊆ Str(J k ), i.e., F consists only of projections and partial projections.

Lemma 7. Let D 0 , D 1 , . . . be as in Lemma 6. Then, for each ≥ 0, D does not contain any partial projection p of arity t with R t ⊆ dom(p).

Proof. We prove Lemma 7 by induction on ≥ 0. For = 0, we have

D 0 = (P \ {p t }) ∪ J k ,
to show that if f ∈ D (n) 0 , n ≥ 1, g 1 , . . . , g n ∈ D (t) and p = f [g 1 , . . . , g n ],
then either p is a projection or a partial projection with R t ⊆ dom(p).

Since g 1 , . . . g n ∈ D , it follows from the induction hypothesis that for every i = 1, . . . , n, either g i is a total projection or g i is a partial projection with R t ⊆ dom(g i ).

Suppose first that g i is a partial projection for some i = 1, . . . , n. Then R t ⊆ dom(g i ) and thus there is a v ∈ R t \ dom(g i ). But then, by the definition of the composition of partial functions p := f [g 1 , . . . , g n ] is not defined on the tuple v, and thus v ∈ dom(p). This shows that R t ⊆ dom(p).

Suppose now that each g 1 , . . . , g n is a total projection of arity t. If v j := (0, . . . , 0, 1

f ∈ J k is a projection, then p = f [g 1 , . . . ,
j-th position , 0, . . . 0) ∈ dom(p) since g 1 ( v j ) = • • • = g n ( v j ) = 0 and (0, . . . , 0) ∈ dom(f ).
Case 2. Let t < n. Then there are j, k ∈ {1, . . . , n}, j = k, such that g j = g k . Without loss of generality, suppose that g j = g k = e t j . In this case, the t-tuple v j defined in [START_REF] Börner | Maximal Partial Clones with no finite basis[END_REF] does not belong to dom(p), since the n-tuple (g 1 ( v j ), . . . , g n ( v j )) contains at least two entries equal to 1, and thus (g 1 ( v j ), . . . , g n ( v j )) ∈ dom(f ). This completes the proof of Lemma 7.

From Lemma 7, it follows that p t ∈ P \ {p t } for every t ≥ 2, and thus we have the following corollary.

Corollary 8. Let X, Y ⊆ {2, 3, . . . }. If X = Y , then {p i : i ∈ X} = {p j : j ∈ Y } .
Let F J k be the set of all partial clones contained in Str(J k ). In this section we show that the strong partial clone pclone(f ) contains a family of strong partial subclones of continuum cardinality.

We will distinguish two cases: Case 1. The following condition is satisfied by the matrix A:

(*) there is α ∈ k such that the m-tuple (α, . . . , α) T is not a column of A.

Case 2. The partial function f is not constant.

We will explain in next subsection why those two cases cover all possibilities.

4.1.

Matrix A satisfies condition (*). For notational ease, we assume in this section that α = 0, i.e., we assume that (0, . . . , 0) T is not a column of A. Denote by I r the identity matrix of size r ≥ 2 on k, i.e., I r is the r × r matrix with symbol 1 on its diagonal and with 0 elsewhere. Furthermore, consider the following infinite set of positive integers

E := {2m r + 1≤i≤r m i-1 : r ≥ 0}.
Remark 11. Notice that, for each r ≥ 1, we have

2m r + m r-1 + • • • + m + 1 = m(2m r-1 + • • • + m + 1) + 1.
Consequently, if t > with t, ∈ E, then t ≥ m + 1.

We will construct an infinite set of partial functions {f : ∈ E} each of which is in pclone(f ), and with the property that for every

∈ E, f ∈ pclone({f k : k = , k ∈ E}).
For ∈ E, let M be the matrix whose first n columns are (a 1 i , . . . , a 1 i , a 2 i , . . . , a 2 i , . . . , a m i , . . . , a m i ) T , i = 1, . . . , n, and the last m× columns are the columns of the identity matrix I m× .

Thus, M is a matrix with m × rows and n + m × columns. Note also that each of the tuples (a i 1 , . . . , a i n ) is repeated times in M .

M =                   a 1 1 . . . . . .
                 
Let ρ be the relation that consists of all columns of M . Note that ρ is a relation of arity m × over k. Now let f be the partial function of arity n + m × whose domain is the set of all rows of the matrix M , and that is given by f ( v) = f (a i 1 , . . . , a i n ), where (a i 1 , . . . , a i n ) are the first n coordinates of the tuple v.

Example 12. Let k = {0, 1, 2}, and let f be the ternary partial function defined by

f   0 0 1 0 1 1 2 0 1   =   1 0 0   .
Here m = n = 3 and for 7 = 2 3 + 1 ∈ E, we get

f 7                      0 0 1 1 0 . . . 0 0 0 0 1 0 1 . . . 0
2 0 1 0 0 . . . 1 0 2 0 1 0 0 . . . 0 1      1 1 . . . 1 0 0 . . . 0 0 . . . 0 0                      . I 21
Note that this matrix has 7 rows starting with (0, 0, 1), 7 rows starting with (0, 1, 1), and 7 rows starting with (2, 0, 1).

Lemma 13. For each ∈ E, we have that f ∈ pclone(f ) and that f is not a partial projection. Proof. Clearly, f ∈ pPolρ . We show that if t = , then there is no matrix whose rows are rows of the matrix M and whose columns belong to the relation ρ t .

For the sake of contradiction, suppose that such a matrix B exists.

Then B has m×t rows and n+m× columns. Let B [i, * ] , i = 1 . . . m×t, be all the rows of B and B [ * ,j] , j = 1, . . . , n + m × , be all its columns.

Case 1. t < .

For i = 1, . . . , m × t, the i-th row of B has the form For i = 0, 1, . . . , s let r i be the number of rows of M that are repeated i × t times in B, with r 0 being the number of rows in M that appear only once in B. Then counting the rows of B we get the equality

B [i, * ] = (
r 0 + t × r 1 + 2t × r 2 + 3t × r 3 + • • • + st × r s = m × t,
i.e., r 0 is a multiple of t. As m× is the total number of rows in M , and since at least one row of M is repeated in B, we have that r 0 < m × .

Furthermore, as t ≥ m × + 1, we conclude that r 0 < m × < t and thus r 0 = 0. Hence, every row of M appears i × t times in B with i ≥ 1. It thus follows that B has no column with one 1 and m × t -1 0's, i.e., every column of B is one of the first n columns of M t .

Let B R be the block matrix consisting of the last m × columns to the right of the matrix B. Let us count the number of 0's in B R . Since the relation ρ t contains no zero tuple, every column of B R contains at least t 1's, i.e., every column in B R contains at most mt -t 0's. Given that there are m × columns in B R , there are at most m × × (m × t -t) 0's in B R . Furthermore, every row of B is a row of M , and thus every row of B R contains exactly one 1 and m × -1 0's.

Hence, there are

m × t × (m × -1) 0's in B R . Since m × t × (m × -1) > m × × (m × t -t),
we conclude that at least one column of B R must contain only 0's, which again yields a contradiction.

Corollary 15. If X ⊆ E and ∈ X, then f ∈ pclone({f t : t ∈ X}).

Proof. Since ∈ X, by Lemma 14 we have that f t ∈ pPolρ for every t ∈ X, and thus pclone({f t : t ∈ X}) ⊆ pPolρ . Since f ∈ pPolρ , we have f ∈ pclone({f t : t ∈ X}). For ∈ E, let f be the ( + 1)-ary partial function whose domain is the set of rows of the matrix

M =     1 1 0 0 . . . 0 1 0 1 0 . . . 0 . . . . . . . . . . . . . . . 1 0 0 0 . . . 1     and let f ( v) := 0 for all v ∈ dom(f ). Then f is a subfunction of f [e +1 1 ]
and so f ∈ C. By Corollary 15 we have that for every

X, Y ⊆ E, X = Y , pclone({f t : t ∈ X}) = pclone({f t : t ∈ Y }).
Hence the strong partial clone C contains each member of the family of strong partial clones {pclone({f t : t ∈ X}) : X ⊆ E}, which is of continuum cardinality.

Remark 18. The construction presented in the above section cannot be used to show that pclone(f ) contains a family of strong partial clones of continuum cardinality if the matrix A contains all constant columns. Indeed in such a case the proof of Lemma 14 fails. The following example illustrates this fact.

Example 19. Let k = {0, 1, 2} and let f be the 4-ary partial function given by

f   0 1 2 0 0 1 2 1 0 1 2 2   =   1 2 0   .
Note that here m = 3 and n = 4. Thus E = {2, 7, 22, . . . }. For = 2 we have

M 2 =       
0 1 2 0 1 0 0 0 0 0 0 1 2 0 0 1 0 0 0 0 0 1 2 1 0 0 1 0 0 0 0 1 2 1 0 0 0 1 0 0 0 1 2 2 0 0 0 0 1 0 0 1 2 2 0 0 0 0 0 1

      
Let ρ 2 be the relation of arity 6 that consists of all columns of M 2 . The matrix below shows that f 3 ∈ pPolρ 2 .

f 3       
0 1 2 0 1 0 0 0 0 0 0 0 0 0 1 2 0 0 1 0 0 0 0 0 0 0 0 1 2 1 0 0 0 1 0 0 0 0 0 0 1 2 1 0 0 0 0 1 0 0 0 0 0 1 2 2 0 0 0 0 0 0 1 0 0 0 1 2 2 0 0 0 0 0 0 0 1 0

       =        1 1 2 2 0 0       
Therefore the construction presented in this subsection is not applicable if we start with an n-ary partial function f whose matrix A does not satisfy the condition (*).

4.2.

The partial function f is non-constant. Notice first that any partial function whose matrix does not satisfy the above condition (*) is included in this case. Indeed, let f be such an n-ary partial function.

Then each of the m-tuples (α, . . . , α) T , α ∈ k, is a column of A. If f was a constant function, then f would be a partial projection function, contradicting f ∈ Str(J k ).

Without loss of generality we may assume that the domain of the n-ary partial function f is repetition free. Indeed if there are 1 ≤ i < j ≤ n such that x i = x j for all x = (x 1 , . . . , x n ) ∈ dom(f ), then we can define the (n -1)-ary partial function g by g(x 1 , . . . , x i , . . . , x j-1 , x j+1 , x n ) := f (x 1 , . . . , x i , . . . , x j-1 , x i , x j+1 . . . , x n ) and use g instead of f . Thus we assume that all n columns of A are distinct. We shall make use of the following definitions and notation.

Let Col(A) be the set of all the columns of A, d f := |dom(f )| and let x i+1 if a = x i and i < |X|,

v f := f (A) ∈ k d f .
x 1 if a = x |X| . Furthermore, for x = (x 1 , . . . , x h ) ∈ [ v f ] h and 1 ≤ i ≤ h we define c i ( x) by c i ( x) := (x 1 , . . . , x i-1 , next [ v f ] (x i ), x i+1 , . . . , x h ).
Since the partial function f is non-constant, the set [ v f ] contains at least two different elements, and so c i ( x) = x for all x ∈ [ v f ] h and all i = 1, . . . , h.

Let α ≥ 0 be the number of columns u in the matrix A that satisfy

[ u] = [ v f ].
Without loss of generality we may assume that all those α columns (if any) are the first columns to the left of the matrix A.

Let ≥ 1. We define the relation ρ of arity d f by

ρ := {c i ( v × f ) | 1 ≤ i ≤ d f } ∪ { u × | u ∈ Col(A)}.
Notice that |ρ | = d f + n. We denote by M the matrix with d f rows and whose d f + n columns are the tuples of ρ written in the following order:

c 1 ( v × f ), . . . , c d f ( v × f ), u × 1 , . . . , u × n ,
where u 1 , . . . , u n are the columns of A written in the same order as they appear in A. We denote by f × the ( d f + n)-ary partial function whose domain is the set of all rows of M and defined by

f × (M ) := v × f .
Notice that for every x := (x 1 , . . . , x d f +n ) ∈ dom(f × ), we have that

x 1 , . . . , x d f +α ∈ [ v f ]. Example 20. Let k = {0, 1, 2}, = 3 and f   0 0 0 1 0 1 0 0 2   =   0 0 1   .
Here v f = (0, 0, 1), [ v f ] = {0, 1}, v ×3 f = (0, 0, 0, 0, 0, 0, 1, 1, 1),

A =   0 0 0 1 0 1 0 0 2
Proof. Let h be the ( d f +n)-ary partial function such that dom(h) := k d f × dom(f ) and h(x 1 , . . . , x d f , y 1 , . . . , y n ) := f (y 1 , y 2 , . . . , y n ). Then f × is the restriction of h on the rows of the matrix M , thus f × ∈ pclone(f ). Now since the set [ v f ] contains at least two different elements, we have

that c i ( v × f ) = v × f for all 1 ≤ i ≤ d f . Morever since f is not a partial projection, we have that u × = v × f for all u ∈ Col(A) and so v × f is not a column of M proving that f × ∈ pPolρ .
Next we want to show that f × ∈ pPolρ m whenever = m. In order to do so, we need the following definition. Let M * be the submatrix of

M whose columns are c 1 ( v × f ), . . . , c d f ( v × f ) together with all columns u × such that [ u × ] = [ v f ].
Thus the matrix M * has exactly d f + α columns. Furthermore, let R be the set of all columns of M * , thus Set v ×m f := (v 1 , . . . , v md f ). Let j 1 = j 2 ∈ {1, . . . , md f } and suppose that the j 1 -th row and the j 2 -th row of M * are identical. Now if v j 1 = v j 2 , then none of the tuples c i ( v ×m f ) (i = 1, . . . , d f , i = j 1 and i = j 2 ) can appear as a column of M * because the j 1 -th (j 2 -th)

R = {ω ∈ ρ | [ω] = [v f ]}.
coordinate of c i ( v ×m f ) is v j 1 (is v j 2 )
. This contradicts the fact that all d f + α columns of the matrix M * are distinct.

Therefore we assume that if j 1 , j 2 ∈ {1, . . . , md f } are such that the j 1th row and the j 2 -th row of M * are identical, then v j 1 = v j 2 and in such a case neither c j 1 ( v ×m f ) nor c j 2 ( v ×m f ) appear as columns in M * . Since at least md f -d f rows of M * are repeated, we conclude that at least md f -d f + 1 tuples of the form c i ( v ×m f ) cannot appear as columns of M * . Consequently, the number of all possible distinct columns in M * is at most

md f + α -(md f -d f + 1) = d f -1 < d f ,

  (k) the set of all n-ary partial functions on k and we set Par(k) := n≥1 Par (n) (k).

Definition 3 .

 3 A partial clone C on k is said to be strong if it is closed under taking subfunctions, i.e., if the condition g ≤ f =⇒ g ∈ C holds for every f ∈ C and g ∈ Par(k). For F ⊆ Par(k), we denote its strong closure by Str(F ), i.e., Str(F ) := {g ∈ Par(k) : g ≤ f for some f ∈ F }.

First

  , we show that (I) minimal strong partial clones (i.e., atoms of L Str(P k ) ) do not exist, and (II) the family {D : D is a partial clone and D ⊆ Str(J k )} is of continuum cardinality. Since every strong partial clone contains Str(J k ), it follows that any strong partial clone contains a family of partial subclones of continuum cardinality. However, partial clones strictly included in Str(J k ) are not strong. Hence, in the last section of the paper, we show that (III) every strong partial clone C = Str(J k ) contains a family of strong partial subclones of continuum cardinality. Some results established in this paper have been presented at the 47 th IEEE International Symposium on Multiple-Valued Logic (Novi Sad, Serbia) and have already appeared in the conference paper [4].

where b 1 1 . . . a m n b m b m

 11 , b 2 , . . . , b m are the values of the partial function f . We there-

  Let f be the restriction of g to the domain {(a 1 , . . . , a n )} and let f (a 1 , . . . , a n ) = b with b ∈ {a 1 , . . . , a n }. Choose c = b and define the (n + 2)-ary partial function h by dom(h) = {(a 1 , . . . , a n , c, b), (a 1 , . . . , a n , b, c)}, and h(a 1 , . . . , a n , c, b) = h(a 1 , . . . , a n , b, c) = b and proceed as in Case 1. Thus we have shown the following result. Theorem 5. Let C = Str(J k ) be a strong partial clone. Then there is a strong partial clone C such that Str(J k ) C C. Hence there are no strong minimal partial clones on k.

Corollary 9 . 4 .

 94 The set F J k is of continuum cardinality on k. Now as every strong partial clone on k contains the set Str(J k ) we immediately have the following result. Theorem 10. Let C be a strong partial clone on k. Then the interval of partial clones [J k , C] is of continuum cardinality. Note that none of the partial clones contained in Str(J k ) is strong. Therefore, Theorem 10 does not provide any insight into the interval of strong partial clones [Str(J k ), C]. Hence, for a strong partial clone C, it is natural to ask what is the size of the interval Clone(C) := {D : D is a strong partial subclone of C}. This is discussed in the next section. Strong partial clones generated by a single partial function Let f be an n-ary partial function not in Str(J k ). Then there are m tuples α 1 , . . . , α m ∈ k n , with α i := (a i 1 , . . . , a i n ), such that f such that (b 1 , . . . , b m ) T is not a column of the matrix on the left hand side of the equation. Denote this m × n matrix by A. In the sequel we assume that the n-ary partial function f , the matrix A and the number m of rows of A are fixed.

Proof.

  It is easy to see that f ∈ pclone(f ). Since f is not a partial projection, none of the first n columns of M is equal to (b 1 , . . . , b 1 , b 2 , . . . , b 2 , . . . , b n , . . . , b n ) T . Moreover, each of the last m × columns of M consists of one 1 and m -1 0's, and thus none of those columns can be equal to (b 1 , . . . , b 1 , b 2 , . . . , b 2 , . . . , b n , . . . , b n ) T . Lemma 14. If , t ∈ E, then f ∈ pPolρ t if and only if = t.

  Now let C be a strong partial clone and suppose that C contains an nary partial function f ∈ Str(J k ) such that |dom(f )| = m. Suppose also that f satisfies condition ( * ). Reconsider the sets E and {f : ∈ E} defined above. By Corollary 15, for every X, Y ⊆ E with X = Y , the two strong partial clones pclone({f t : t ∈ X}) and pclone({f t : t ∈ Y }) are distinct and both are contained in pclone(f ) ⊆ C. Therefore, we have the following corollary. Corollary 16. Let C be a strong partial clone on k and suppose that C contains a partial function f ∈ Str(J k ) whose matrix A satisfies ( * ). Then the set of strong partial clones contained in C is of continuum cardinality. Example 17. Let C be a strong partial clone on k and suppose that C contains a unary partial function g such that 1 ∈ dom(g) and g(1) = 0. Let f := g| {1} , i.e., f is given by dom(f ) = {1} and f (1) = 0. Since C is a strong partial clone, we have f ∈ C. Consider the set E defined above with m = 1, i.e., E = {2, 3, . . . }.

  For x := (x 1 , . . . , x h ) ∈ k h and ≥ 1, we define x × ∈ k h by x × = (x 1 , . . . , x 1 times , x 2 , . . . , x 2 times , . . . , x h , . . . , x h times ), and let [ x] := {x 1 , . . . , x h }. For a set X ⊆ k with X = {x 1 < • • • < x |X| } and a ∈ X we define next X (a) ∈ X by next X (a) :=

  Denote by F the ( d f + α)-ary partial function with domain the set of all rows of M * and defined by F (M * ) := v × f . Now by construction of f × , we have that the value of f × (x) is uniquely determined by the first d f entries of x for every row x of the matrix M . Therefore the value of f × is determined by the first d f columns of M and thus the partial function F is well defined. Summing up, we have the following result.Lemma 23. For all , m ≥ 2, if F ∈ pPolR m then f × ∈ pPolρ m .Proof. Suppose that f × ∈ pPolρ m . So there is a matrix P whose d f + n columns are all from ρ m and whose rows are from dom(f × ), and such that f × (P ) / ∈ ρ m . Let P be the submatrix of P consisting of the first d f + α columns of P .Since the firstd f + α entries of every tuple x ∈ dom(f × ) belong to [v f ],we have that all columns of P belong to R m , and since F (P ) = f × (P ) / ∈ ρ m we have that F ∈ pPolR m .Lemma 24. Let , m ≥ 2 with = m. Then f × ∈ pPolρ m .Proof. We distinguish two cases. Case 1. m < . Let M be a matrix whose rows are in dom(f × ) and whose columns are in ρ m . Thus M has md f rows all of which are rows of M and M has d f + n columns all of which are columns of M m . As m < the matrix M has less rows than M and so by Remark 21, f × restricted to the rows of M is a partial projection function. Thus f × (M ) ∈ ρ m proving that f × ∈ pPolρ m . Case 2. < m. We argue the contrapositive. Suppose that f × ∈ pPolρ m . Then there is a matrix M whose rows are in dom(f × ) and whose columns are in ρ m and such that f × (M ) ∈ ρ m . If some row of M is missing in M , then by Remark 21, f × restricted to the rows of M is a partial projection function, but this contradicts our assumption that f × (M ) ∈ ρ m . Therefore we assume that all rows of M appear in M and consequently the matrix M contains md f -d f repeated rows. Also, as all columns of M are distinct, we have that all d f + n columns of M are distinct. Let M * be the submatrix of M consisting of the first d f + α columns of M . Since f × (M ) ∈ ρ m , we have, by Lemma 23, that F (M * ) ∈ R m . Since all columns of M are distinct, we have that all d f + α columns of M * are distinct. Furthermore, since the matrix M contains md f -d f repeated rows we have that M * contains at least md f -d f repeated rows.

  and thus the matrix M * cannot have all its d f columns distinct. This contradiction shows that F ∈ pPolR m , and by Lemma 23f × ∈ pPolρ m .The proof of the following result is similar to that of Corollary 16.Corollary 25. Let C be a strong partial clone on k and suppose that C contains a partial function f ∈ Str(J k ) that is not a constant function. Then the set of strong partial clones contained in C is of continuum cardinality.

  g n ] is a total projection and we are done. So suppose that f is a partial projection. Since f ∈ D 0 , we have f = p

n with n = t.

Case 1. Let n < t. Then there is a j ∈ {1, . . . , t} such e t j ∈ {g 1 , . . . , g n }. Therefore

[START_REF] Börner | Maximal Partial Clones with no finite basis[END_REF]

  Now if M [i, * ] appears repeated exactly r times in B, 2 ≤ r ≤ m × t, then B [ * ,n+i] will contain r 1's. Given that this column belongs to ρ t and given that each row of A occurs t times in M t , we conclude that r is a multiple of t.

	Case 2. < t.
	As , t ∈ E, we must have t ≥ m × + 1, and
	M [i, * ] = (a s 1 , . . . , a s n , 0, . . . , 0,	1	, 0, . . . , 0)
		i-th position

a s 1 , . . . , a s n , x 1 , . . . , x m× ) where s ∈ {1, . . . , m} and exactly one of x 1 , . . . , x m× is equal to 1 and all other are 0. Given that B consists of m × t such rows and given that mt < m × , there is a j ∈ {n + 1, . . . , n + m × } such that the column B [ * ,j] consists of only 0's. Hence, it does not belong to ρ t , which yields the desired contradiction. where 1 ≤ s ≤ m. Therefore, if for some i ∈ {1, . . . , m × }, the i-th row M [i, * ] of M does not appear in B, then B [ * ,n+i] will contain only 0's, and it will not belong to the relation ρ t . Thus every row of M must appear in the matrix B. Since < t, some of the rows of M are repeated in B.