On an ODE-PDE coupling model of the mitochondrial swelling process

Mitochondrial swelling has huge impact to multicellular organisms since it triggers apoptosis, the programmed cell death. In this paper we present a new mathematical model of this phenomenon. As a novelty it includes spatial eﬀects, which are of great importance for the in vivo process. Our model considers three mitochondrial subpopulations varying in the degree of swelling. The evolution of these groups is dependent on the present calcium concentration and is described by a system of ODEs, whereas the calcium propagation is modeled by a reaction-diﬀusion equation taking into account spatial eﬀects. We analyze the derived model with respect to existence and long-time behavior of solutions and obtain a complete mathematical classiﬁcation of the swelling process.


Introduction
Biological background Mitochondria are often termed the cell's powerhouse due to their main function as energy supplier for almost all eukaryotic cells [1]. However, these double-membrane enclosed organelles also play a decisive role in cell death by their ability to trigger apoptosis. One of the key factors in this process is the permeabilization of the inner mitochondrial membrane [13], resulting in the swelling of the mitochondrial matrix. Mitochondrial permeability transition is effectuated by the opening of a pore in the inner membrane, which happens under pathological conditions like high Ca 2+ concentrations [14]. The increased permeability leads to an osmotically driven influx of solutes and water into the mitochondrial matrix, which in turn causes swelling [8], [14]. This process culminates in the rupture of the outer membrane [20]. Outer membrane rupture denotes a critical event, since now apoptosis is irreversibly triggered by the release of several proapoptotic factors from the intermembrane space [15]. Intact mitochondria store calcium in their matrix. If swelling is induced, this stored calcium is additionally released [15] and the remaining mitochondria are confronted with an even higher calcium load, leading to an acceleration of the process.
Experimental procedure The swelling of mitochondria is measured on the basis of light scattering. While intact mitochondria show high light scattering values, the more mitochondria are swollen the less light is deflected. Hence the volume increase is indirectly displayed by a decreasing optical density. This relation is shown to be linear [5], [11], [21]. The experiments were performed at the Helmholtz Zentrum München, Institute of Molecular Toxicology and Pharmacology using rat liver mitochondria and Ca 2+ as swelling inducer.
Existing models Although the process of mitochondrial swelling induced by calcium is known for more than 30 years, mathematical modeling has only started recently. At this there are two conceptually different approaches: The microscale models focussing on a detailed description of all biochemical processes in single mitochondria [22], [25], and the macroscale models which directly aim to represent the swelling of a whole population [5], [11], [17]. In our previous publication [11], we presented a mathematical model that is for the first time capable to describe the whole time progress of mitochondrial swelling. The model is based on the observation that mitochondria vary within subpopulation concerning their sensitivity for swelling induction as it was described in [26]. It consists of one ODE for the fraction of swollen mitochondria and a delay equation to determine the corresponding population volume. This model is in great accordance with the experimentally obtained data and shows consistent parameter values with increasing Ca 2+ concentrations. In [10], it has been shown (see Fig2.1, page 15) the dependence of swelling curves at different calcium concentrations and how the swelling curves can be modeled with the simple mathematical model. The most important parameters of the model among others is the feedback parameter and the parameter depicting the average swelling time of mitochondria. The analysis of the simple model also yielded that the experimental swelling curves can only be described in an accurate way by including the positive feedback, otherwise one could not receive the correct shape of the swelling curves.
Furthermore, different swelling inducers and mitochondria from different organs can be classified by comparing the corresponding optimal parameter values.
Spatial effects Experiments revealed the necessity to take into account spatial effects. In fact, models taking into account of spatial effects are discussed in [2] and [3].
The same amount of Ca 2+ is added in different concentrations, i.e. mitochondria are exposed to identical calcium quantities with varying levels of concentration profiles. This leads to different shapes of the corresponding swelling curves, which only trace back to the different calcium distributions. Obviously this implies the influence of spatial effects. In particular the dependence on local processes gets important when thinking of mitochondrial swelling taking place in vivo. There are two mechanisms that lead to intracellular Ca 2+ increase (see e.g. [23], [24]): Internal release from the endoplasmic reticulum or exter-nal calcium influx from the extracellular milieu. Both calcium sources are highly localized. Furthermore mitochondria within cells are not distributed randomly but reside in three main regions [18], which enforces the influence of spatial effects.

The mitochondria model
In this work we introduce a new mathematical model that takes into account the above mentioned spatial effects. This results in a coupled ODE-PDE system.

Description
In accordance with our theoretical [11] and experimental [26] findings, we assume that three subpopulations of mitochondria with different corresponding volumes exist. Here N 1 (x, t) describes the density of intact, unswollen mitochondria, N 2 (x, t) contains mitochondria that are in the swelling process but not completely swollen and N 3 (x, t) denotes the density of completely swollen mitochondria. The Ca 2+ concentration is denoted by u(x, t). The transition of intact mitochondria over swelling to completely swollen ones proceeds in dependence on the local calcium concentration. At this we can assume that mitochondria do not move in any direction and hence the spatial effects are only introduced by the calcium evolution. The evolution of the mitochondrial subpopulations is modeled by a system of ODEs, that depends on the space variable x in terms of a parameter. We analyze the swelling of mitochondria on a bounded domain Ω ⊂ R n with n = 2, 3. This domain could either be the test tube or the whole cell. The initial calcium concentration u(x, 0) describes the added amount of Ca 2+ to induce the swelling process.
This leads to the following coupled ODE-PDE system determined by the non-negative model functions f and g: with diffusion constant d 1 > 0 and feedback parameter d 2 > 0. The boundary conditions are given by ∂ ν u = 0 on ∂Ω (5) and we have initial values Model function f The process of mitochondrial permeability transition is dependent on the calcium concentration. If the local concentration of Ca 2+ is sufficiently high, the pores open and mitochondrial swelling is initiated. This incident is mathematically described by the transition of mitochondria from N 1 to N 2 . The corresponding transition function f (u) is zero up to a certain threshold C − displaying the concentration which is needed to start the whole process. Whenever this Ca 2+ threshold is reached, the local transition at this point from N 1 to N 3 over N 2 is inevitably triggered. It is written e.g. in [21] that this process is calcium-dependent with higher concentrations leading to faster pore opening. Hence the function f (u) is increasing in u. The transfer from unswollen to swelling mitochondria is related to pore opening, hence we also postulate that there is some saturation rate f * displaying the maximal transition rate. This is biologically explained by a bounded speed of pore opening with increasing calcium concentrations.

Remark
The initiation threshold of f is crucial for the whole swelling procedure. Dependent on the amount and location of added calcium, it can happen that in the beginning the local concentration was enough to induce swelling in this region, but after some time due to diffusion the threshold C − is not reached anymore. Thus we only have partial swelling and after the whole process there are still intact mitochondria left. Nevertheless, there are no mitochondria in the intermediate state N 2 .
Model function g The change of the population N 2 consists of mitochondria entering the swelling process (coming from N 1 ) and mitochondria getting completely swollen (leaving to N 3 ). The transition from N 2 to N 3 is modeled by the transition function g(u). In contrast to the function f , here we have no initiation threshold and this transition can not be avoided. This property is based on the biological mechanism. The permeabilization of the inner membrane due to pore opening leads to water influx and hence unstoppable swelling of the mitochondrial matrix. This process itself is independent of the present calcium concentration. Due to a limited pore size, this effect also has its restriction and thus we have saturation at level g * . However, biologically it is not clear if there are other influences of calcium to this second transition, e.g. by the opening of additional pores. To include such possibilities, we allow for general increasing g with saturation at level g * .
In [26], they proved that there are different mitochondrial subpopulations and gave the shapes of the functions f and g as described in the section "Numerical Simulations" discussed with biologists.
The third population N 3 of completely swollen mitochondria grow continuously due to the unstoppable transition from N 2 to N 3 . All mitochondria that started to swell will be completely swollen in the end. Calcium evolution The model consists of spatial developments in terms of diffusing calcium. In addition to the diffusion term, the equation for the calcium concentration contains a production term dependent on N 2 , which is justified by the following: In our earlier ODE approach [11] we showed that it is essential to include a positive feedback mechanism. It is not possible to display the correct swelling behavior without the positive feedback, even with the much more simple model. In the biological community there is no doubt about the existence of the positive feedback. This accelerating effect is induced by stored calcium inside the mitochondria, which is additionally released once the mitochondrion gets completely swollen. Due to a fixed amount of stored Ca 2+ , we assume that the additionally released calcium amount is proportional to the newly completely swollen mitochondria, i.e., those mitochondria leaving N 2 and entering N 3 . Here the feedback parameter d 2 describes the amount of stored calcium.
In this paper we are interested in the well-posedness and long-time behavior of solutions. The unique existence of the global solution is obtained by the contraction mapping principle. Furthermore we present a classification of limiting profiles. Depending on an a-priori given threshold, we show two possible scenarios, i.e., partial and complete swelling.

Well-posedness
The coupled ODE-PDE model (1) -(5) describing the mitochondrial swelling process will now be analyzed mathematically. At first we want to show the global existence and uniqueness of the solution (u, N 1 , N 2 , N 3 ) on the phase space L 2 (Ω). For that purpose we introduce some assumptions to the model functions f and g.

Condition 1
For the model functions f : R → R and g : R → R it holds: (i) Non-negativity and Boundedness: (ii) Lipschitz continuity: One remarkable characteristic of the model is the following conservation law: i.e. the total population N := N 1 (x, t) + N 2 (x, t) + N 3 (x, t) does not change and is given by the sum of the initial data. In fact, adding three equations (2) + (3) + (4), we obtain ∂ t N = 0. The following theorem yields the desired result of well-posedness.

Theorem 1
Let Ω be a bounded domain in R n and let Condition 1 be satisfied. Then for all initial data u 0 ∈ L 2 (Ω), N i,0 ∈ L ∞ (Ω) (i = 1, 2, 3), the system (1) -(5) possesses a unique global solution (u, Proof We first note that by (6) the essential unknown functions can be taken as (u, N 1 , N 2 ). Let X T := C([0, T ]; L 2 (Ω)) and define the mapping Here for a given u ∈ X T , N u = (N u 1 , N u 2 ) denotes the solution of the ODE problem : andû denotes the solution of the pure PDE problem : Obviously, by Condition 1 F u is Lipschitz continuous from Y = L ∞ (Ω) × L ∞ (Ω) into itself, the Picard-Lindelöf theorem assures the existence of the unique global solution is Lipschitz continuous from L 2 (Ω into itself by Condition 1, the standard argument shows that (8) has the unique solutionû ∈ C([0, ∞); loc ((0, ∞); L 2 (Ω)) ( see e.g. [4], [9] or in a more abstract way in [6], [19] ). In order to show that B becomes a contraction mapping in X T for a sufficiently small T ∈ (0, 1], we are going to establish a priori estimates for the difference of two solutions. We put Then by virtue of the boundedness of f and g, we obtain Step 1 Multiplying (9) by δN 1 and using the Lipschitz continuity and positivity of f , we get d dt Hence by (12) and (13), we get Step 2 Multiplying (10) by δN 2 and using the Lipschitz continuity of f, g, the positivity of g, boundedness of f , (13) and (14), we now have Then by (12) and (15), we get Step 3 Multiplying (11) by δu and using the Lipschitz continuity and boundedness of g, we have Hence by (16) and (12), we obtain for any t ∈ (0, S] and S ∈ (0, T ] Consequently we find Thus there exists a sufficiently small T 0 depending on T and other parameters but not on the initial data such that B possesses a unique fixed point u ∈ X T 0 . In other words, gives a solution of the system (1) - (5). The uniqueness of (N u 1 , N u 2 ) follows from (15) and (16) directly. Since T 0 > 0 does not depend on the choice of initial data, it is easy to see that this local solution can be continued up to [0, T ] for any T .

2
The model variables have a biological meaning and hence it is important to show that they are non-negative. This is done in the following.

Proposition 2
Let all assumptions of Theorem 1 hold and in addition assume that Then the solution (u, N 1 , N 2 , N 3 ) preserves non-negativity. Furthermore N 1 , N 2 and N 3 are uniformly bounded in Ω × [0, ∞).

Proof 1.) Non-negativity
Multiply the equations by the negative part of solutions (u . Then by using the Lipschitz continuity of f, g and Gronwall's inequality, we can deduce From the conservation law (6) and the proved non-negativity of the ODE parts N 1 , N 2 and N 3 it follows immediately

Asymptotic behavior of solutions
Now the longtime behavior of the solution (u, N 1 , N 2 , N 3 ) is studied. This behavior is highly dependent on the special structure of the model functions f and g.

Proposition 3
Let all assumptions of Theorem 1 and Proposition 2 hold and in addition assume u 0 ≡ 0. Then the unique solution u is strictly positive for t > 0 and becomes bounded below by a strictly positive constant:

Proof
The solution u fulfills the PDE problem For the proof we introduce a subsolution u of (18) by Since d 2 g(u)N 2 ≥ 0, it follows from the comparison principle that Comparing u with the subsolution u ≡ 0 of (19) and applying the strong parabolic maximum principle, we obtain the strict bound u(x, t) > 0. Furthermore, we note that Hop's maximum principle together with the fact u(x, t) > 0 and the homogeneous Neumann boundary condition assures the strict positivity of u(x, t) > 0 on the boundary. Thus we find that Then u ≡ satisfies Again applying the comparison principle, we obtain This result is now used to obtain information about the type of convergence as time goes to infinity. For that we need additional assumptions on the functions f and g.

Condition 2
Let the model functions f : R → R and g : R → R fulfill Condition 1. In addition we assume that there exist constants C − > 0, m 1 > 0, m 2 > 0, δ 0 > 0 and 0 > 0 such that the following assertions hold: (i) Starting threshold: (iii) Lower bounds: In order to show several convergence results we need the following proposition, which can be easily derived from the standard Gronwall's Lemma.

Proposition 4
Let y(t) and a(t) be non-negative functions with y ∈ C 1 ([t 0 , t 1 )) and a ∈ C([t 0 , t 1 )) for 0 ≤ t 0 < t 1 ≤ ∞. Suppose that y satisfies for some γ 0 > 0 d dt Then the following estimates hold true.
Multiplying the differential inequality by e γ 0 t , we easily get Hence with the aid of simple calculations, we can deduce statements above.
The next theorem gives information about the strong convergence of the solution.

Theorem 5
Let Condition 2 and assumptions in Proposition 2 be satisfied. Then we have the following strong convergence results:

Proof
(1) From the model equation (2) and the non-negativity result it holds in the point-wise Hence the sequence is non-increasing and bounded below by 0, whence follows the convergence Furthermore, by (17) we find Then by virtue of the Lebesgue dominated convergence theorem, we conclude that N 1 (·, t) converges to N ∞ 1 (·) strongly in L 1 (Ω) as t → ∞. Thus to deduce (21), it suffices to use the relation (2) As for N 3 (x, t), the model equation (4) gives Since N 3 (x, t) is bounded above by N L ∞ (Ω) , the monotonicity yields the almost everywhere convergence Moreover |N 3 (x, t)| is also dominated by N L ∞ (Ω) , then by the same arguments for N 1 above we can deduce (23).
(3) Combining (6) with (25) and (26), we easily obtain the almost everywhere convergence for for a.e. x ∈ Ω . and the convergence in L p (Ω) follows immediately as before.
(4) Here we are going to show that N ∞ 2 (x) ≡ 0. To this end, we first note that the integration of (4) on (0, t) gives Then by (28), Proposition 3 and (iii), (iv) of Condition 2, it holds Therefore there exists a sequence {t k } k∈N with t k → ∞ such that whence follows N ∞ 2 (x) ≡ 0 for a.e. x ∈ Ω . (5) In order to show the convergence properties of u, we use the following orthogonal decomposition where ϕ 1 (x) ≡ C ϕ = |Ω| −1/2 is the eigenfunction for the first eigenvalue λ 1 = 0 of −∆ with the homogeneous Neumann boundary condition and ϕ ⊥ (x, t) is orthogonal to ϕ 1 (x), i.e., Multiplying (1) by ϕ 1 (x) and using the integration by parts, we get Hence the function a 1 (t) is non-decreasing in t. Furthermore, substituting the relation g(u(x, t))N 2 (x, t) = ∂ t N 3 (x, t) into (32), we obtain Thus we find In order to show that u converges to a constant function, it is thus sufficient to show that ϕ ⊥ (t, x) → 0 as t → ∞ for a.e. x ∈ Ω. Here by Wirtinger's inequality (see [7]), we get Furthermore, for any v ∈ H ⊥ , using the fact that v satisfies the homogeneous Neumann boundary condition and Wirtinger's inequality again, we get Thus we obtain Substitution of the decomposition (30) into the PDE (1) leads to Multiplying (36) by ϕ ⊥ , we get Then by (35) and Hölder's inequality, we obtain Here we note (29) implies Then applying (ii) of Proposition 4 with which implies (24) We can obtain further estimates:

Proposition 6
Under the assumptions of Theorem 5 the following additional facts hold:

Classification of partial and complete swelling
The mitochondrial swelling process and its extent is dependent on the local calcium dose. If the initial concentration u 0 stays below the initiation threshold C − at any point x ∈ Ω and N 2,0 = 0, then no swelling will happen and we have N i (x, t) ≡ N i,0 (x) ∀x ∈ Ω, i = 1, 2, 3 for all t > 0. Another possible scenario is the so-called "partial swelling". This effect of partial swelling occurs in the experiments and can also be seen in the simulations when the initial calcium concentration lies above C − at a small region but due to diffusion it does not stay above this threshold for the whole time. This leads to the case where there exists a finite time But if the initial calcium distribution together with the influence of the positive feedback is sufficiently high, then "complete swelling" occurs which means N 1 (x, t) → 0 and N 3 (x, t) → N (x) for all x ∈ Ω as t → ∞.
As it was shown before, for both cases it holds that N 2 (x, t) → 0 as t → ∞.

Condition 3
Let the assumption of Conditions 1 and 2 hold. In addition we assume more regularity of the initial data: A crucial point to distinguish between partial and complete swelling is to check if f (u) stays positive for all times. For that it is necessary to have uniform convergence of u(x, t) to u ∞ ≡ a ∞ 1 C ϕ . Up to now we only have the strong convergence in L 2 (Ω). So our aim now is to show the uniform convergence, which turns out to be an extensive task.

Theorem 7
Let N ≤ 3 and Condition 3 be satisfied, then the following additional statements hold:
which assures that the second term of the right-hand side of (50) is bounded by some constant C 4 for all t ∈ (0, ∞). Furthermore, since by Theorem 1 the solution u fulfills √ t ∆u ∈ L 2 (0, T ; L 2 (Ω)), for each δ > 0, there exist a constant C 5 and δ ∈ (0, δ) such that Moreover (43) assures that the first term of the right-hand side of (50) is bounded by some constant C 6 for all t ∈ (δ, ∞). Thus in order to deduce (48), it suffices to apply Proposition 4 to (50) with Then in order to show the uniform convergence of u(t) to u ∞ , it suffices to verify (51). For this purpose, we begin with the estimate for ∇N 1 (t).
Application of the gradient to the model equation (2) leads to Here we note that since f is Lipschitz and u ∈ H 1 (Ω), f (u)∇ exist a.e. x ∈ Ω and belongs to H 1 (Ω) ( see e.g. [12] and [16] ). Then multiplying (52) by ∇N 1 (t), we get Here we recall In view of Condition 2, we decompose Ω into 3 parts : Then (i) of Condition 2 implies that f (u(x, t)) = 0 in Ω 1 (t). According to the behavior of f (u(x, t)) in Ω 2 (t), we have to distinguish between the two cases Case (I) We first take a look at the case u ∞ ≤ C − .
On Ω 2 (t), in view of (ii) of Condition 2 and (54), we have Substituting this estimate into (53), we get The first integral can be further estimated as follows : Here we used the Sobolev embedding theorem, the elliptic estimate in L 2 (Ω) and Wirtinger's inequality ( C H 1 is the embedding constant for H 1 (Ω) ⊂ L 4 (Ω) and C 7 is a constant depending on C H 1 , C W ).
For the second integral over Ω 3 (t), recalling that f (u(x, t)) ≥ f (C − + δ 0 ) > 0 in Ω 3 (t) by (iii) of Condition 2, we artificially insert and get Substitution of these findings into (55) leads to The last term can be omitted and by Young's inequality, there exists a constant C 8 such that Integration of this over (δ, t) yields Here we note that (53) gives whence follows Thus, in view of (42) and (43) we deduce from (57) and (56) that This result is now used to show that ∇N 2 (t) also stays bounded in L 2 (Ω). For that we apply the gradient to model equation (3) and get Multiplication by ∇N 2 and integration over Ω leads to Then (ii), (iii) of Condition 1 yields Hence by Young's inequality and (58), for any η > 0, there exists a constant C 9,η such that Then (59) with η = 1 and Gronwall's inequality give where t * > 0 is the number given in Proposition 3, from which we know that for all x ∈ Ω it holds u(x, t) ≥ > 0 ∀t ≥ t * and consequently (iii) and (iv) of Condition 2 imply g(u(x, t)) ≥ g( ) > 0 ∀t ≥ t * with = min( , 0 ) as defined earlier. Substituting this into the last term of (59) with η = g( ), we get Then in view of (42) and (43), we can apply (i) of Proposition 4 with Thus together with (60), we achieved (51).
Case II We now take a look at the case u ∞ > C − . Let x ∈ Ω 2 (t), i.e., C − ≤ u(x, t) ≤ C − + δ 0 , then recalling (ii) of Condition 2 : and integrating this over (C − , u(x, t)) with respect to s, we get Due to the condition u ∞ > C − and the monotonicity of a 1 (t) with (54) it follows: For all α with 0 < α < (u ∞ − C − ), there exists T α > 0 such that Substituting this into (61), we obtain with Young's inequality Moreover since this estimate is still valid for the case where u( Since f (u(x, t)) = 0 = f (u(x, t)) in Ω 1 (t) for t ≥ T α , substituting these estimates for f (u(x, t)) into (53), we have Hence by Young's inequality, Wirtinger's inequality and the Sobolev space embedding H 2 (Ω) → L ∞ (Ω), there exists a constant C 10 such that Integrating this over (T α , t) and applying Gronwall's inequality, in view of (42) and (43), we find that Here by the same reasoning as for (57), we get sup t ∈ [0,Tα] ∇N 1 (t) L 2 (Ω) < ∞ . Thus we obtain sup The boundedness of ∇N 2 (t) in L 2 (Ω) can be shown in exactly the same way as for case 1) u ∞ ≤ C − , since for the estimation of ∇N 2 (t) L 2 (Ω) we did not use the relation between u ∞ and C − . This yields sup Thus (51) is verified and the proof is completed. 2 Under the assumptions of Condition 3, the long-time behavior can be further characterized as follows:

Theorem 8
Let Condition 3 be satisfied, then the following asymptotic behavior of solution hold.

1.) Partial swelling
Let u ∞ < C − , then there exists a finite time T p > 0 such that and we have the following exponential convergence rates for all x ∈ Ω and t ≥ T p : These rates are dependent on the model function g and the earlier defined parameter = min( , 0 ) in correspondence to Condition 2 and Proposition 3 and we have

2.) Complete swelling
Let u ∞ > C − , then there exists some T c > 0 such that for all x ∈ Ω and all t ≥ T c the following exponential convergence rates hold.

Proof
1.) Partial swelling We start with the partial swelling case. By the uniform convergence of u(x, t) to u ∞ < C − assures the existence of a time T p > 0 such that For equation (3), by (62) and the definition of , it holds for all t ≥ T p which gives i.e., we have exponential convergence of N 2 (x, t) to 0. From the conservation law we know which together with (63) yields Now, in view of (33) and (64), we can easily see that As for the convergence of ϕ ⊥ (t) 2 L 2 (Ω) , we apply (ii) of Proposition 4 with and t 0 = t 2 to (38). Then we obtain by (63) which implies As for ∇ϕ ⊥ (t) 2 L 2 (Ω) , from (44) we now get Thus repeating the same argument as for ϕ ⊥ (t) 2 L 2 (Ω) , we can obtain the convergence

2.) Complete swelling
In the complete swelling case where u ∞ > C − , the uniform convergence of u(x, t) to u ∞ implies that for any β ∈ (0, u ∞ − C − ), there exists a time In order to derive the strict positivity of f (u(x, t)), we have to distinguish between two cases in accordance with Condition 2: • u(x, t) ≥ C − + δ 0 , where δ 0 > 0 denotes the constant from Condition 2, i.e., β ≥ δ 0 . In this case, it follows from (iii) of Condition 2 that • u(x, t) ∈ (C − , C − + δ 0 ), which means we are in the case (ii) of Condition 2, where we have the relation f (u(x, t)) ≥ m 1 β > 0, whence follows In summary we conclude and in addition by (iii) and (iv) of Condition 2 Substituting relation (65) into the model equation (2), we obtain the exponential decay of N 1 (x, t) to 0: Furthermore the second model equation (3) can be estimated by means of the previous result, (ii) of Condition 1 and (66) as follows.
Then by (67), we easily get Hence integrating (68) over (T c , t), we obtain the exponential decay of N 2 (x, t) : In analogy to the previous case, the conservation law together with the exponential decay obtained above implies The exponential convergence for a 1 (t), ϕ ⊥ (t) 2 L 2 (Ω) and ∇ϕ ⊥ (t) 2 L 2 (Ω) can be derived from the same reasoning as for the partial swelling case with g( ) replaced by η. 2

Numerical simulation
The presented model will now be analyzed numerically. Here we consider the in vitro case, i.e. the domain Ω describes a test tube containing purified mitochondria. We are interested in the effects of adding a certain calcium amount to intact mitochondria. The time development of the calcium concentration u and the mitochondrial subpopulations N 1 (intact), N 2 (swelling) and N 3 (completely swollen) is simulated using MATLAB.
Based on biological observations mentioned in the beginning we assume a sigmoidal shape of the model functions f and g determined in the following way: The model parameters we used for the simulations are noted in Table 1.
The calcium concentrations C − and C + are adapted from [10], the other parameters are chosen exemplary. The step size of space discretization is 1/40 and the number of grid points = 40 2 .

Numerical approximation
The PDE describing the calcium diffusion process is discretized with respect to space by means of the standard finite difference approach. Here the Laplacian is approximated by use of the five point star. Doing so, the PDE is transferred into an ODE and we end up with an ODE system on the discrete domain. Due to the low numerical complexity of the model, this can be easily achieved by using the explicit Euler method. The homogeneous Neumann boundary condition is realized by introducing phantom points in order to calculate the normal derivative at the boundary.

Initial values
As we pointed out earlier, in the beginning all mitochondria are intact and with that neither in the swelling process nor completely swollen, i.e. N 1,0 (x) ≡ 1, N 2,0 (x) ≡ 0, N 3,0 (x) ≡ 0 .
For the calcium concentration it is not so clear how to determine the initial state. The initial value u 0 (x) defines the distribution of the added Ca 2+ amount. At this the rate of diffusion progression as well as the dosage location is of great importance. Therefore one can imagine different possible initial states for a fixed total amount C tot of added calcium. For the simulations we used C tot = 30 · (N + 1) 2 dependent on the space discretization and the initial calcium distribution is assumed to be the normal distribution restricted to the sector [−1, 3] projected onto the discrete domain. For the calcium concentration, it is not clear how to determine initial state. Therefore the initial calcium concentration is determined by a sector of the standard normal distribution (see [10], pages 89,90]. ) The simulations show that this calcium amount is high enough to induce complete swelling.

Results: Complete Swelling
The columns in Figure 3 show the evolution of the model variables u, N 1 , N 2 and N 3 .

Remark
One has to be aware of different time scales. The calcium diffusion occurs slower than the development of the mitochondrial populations and the constant state is reached later (t = 4100 for the calcium evolution, t = 1900 for the mitochondrial populations).
Based on the initial calcium distribution swelling does not start on the whole domain immediately, but only on the regions where the concentrations exceeds the swelling threshold C − . One remarkable result is the clearly visible spreading calcium wave, which particularly becomes obvious in the evolution of N 2 . If we compare the dynamics with those of simple diffusion without any feedback, the resulting calcium evolution induced by mitochondrial swelling is indeed completely different. In accordance with the analytical results, in the end all mitochondria are completely swollen and the calcium concentration is constant on the whole domain.

Results: Partial swelling
Now we consider a much lower total amount of calcium, C tot = 2.1 · (N + 1) 2 added in two varying initial distributions. The simulations show that this small change in the degree of localization decides between complete and partial swelling. In Figures 4 and 5 we only show the evolution of u and N 1 in order to get a qualitative description of the difference between complete and partial swelling , by which it is shown how a slightly modified initial condition can lead to partial swelling (on the right side) or complete swelling (left).

Conclusion
The developed mathematical model is in total accordance with the biologically expected results. It provides a deep understanding of the underlying mechanism with focus on the spatial development. This is of great importance for the understanding of the processes taking place in vivo.
The mathematical analysis yields interesting results and in particular we were able to obtain a complete classification of the mitochondrial swelling process. The robustness and validity of the derived model were shown.
This model can also be adapted to the swelling process taking place in cells instead of test tubes by changing the boundary condition. A cell is not a closed system anymore and hence there is a calcium flux over the boundary leading to Robin boundary conditions.