A parabolic-hyperbolic system modeling the growth of a tumor
AA PARABOLIC-HYPERBOLIC SYSTEM MODELING THE GROWTH OF ATUMOR
RUI LI AND BEI HU
Abstract.
In this paper, we consider a model with tumor microenvironment involving nutrientdensity, extracellular matrix and matrix degrading enzymes, which satisfy a coupled system of PDEswith a free boundary. For this coupled parabolic-hyperbolic free boundary problem, we prove thatthere is a unique radially symmetric solution globally in time. The stationary problem involvesa ODE system which is transformed into a singular integro-differential equation. We establish awell-posed theorem for such general types of equations by the shooting method; the theorem isthen applied to our problem for the existence of a stationary solution. In addition, for this highlynonlinear problem, we also prove the uniqueness of the stationary solution, which is a nontrivialresult. In addition, numerical simulations indicate that the stationary solution is likely locallyasymptotically stable for reasonable range of parameters. Introduction
It is estimated that there are 8.2 million cancer-related deaths worldwide every year. Tumormalignancy and metastatic progression are the primary cause, which leads to 90 percent of deathsfrom cancer. Many recent cancer-related studies have pointed out that the remodeling of collagenfibers in the extracellular matrix (ECM) of the tumor microenvironment facilitates the migrationof cancer cells during metastasis, since such modifications of ECM collagen fibers result in changesof ECM physical and biomachanical properties that affect cancer cell migration through the ECM[27]. The ECM is defined as the diverse collection of proteins and sugars that surrounds cells inall solid tissues. This tissue compartment provides structural support by maintaining an insolublescaffold, and this in turn defines the characteristic shape and dimensions of organs and complextissues [22]. Actually, various types of fibrous proteins are present in the ECM including collagens,elastins and laminis; among these, collagen is the most abound-ant and the main structural proteinin the ECM [3]. In general, the ECM degradation caused by enzyme matrix is a key procedurefor the ECM remodeling. In this paper, we try to use the matrix degrading enzymes (MDE) todescribe the degrading process.Over the last few decades, mathematical modeling has played a vital role in testing hypotheses,simulating the dynamics of complex systems and understanding the mechanistic underpinningsof dynamical systems. In particular, an increasing number of mathematical models describingsolid tumor growth have been studied and developed; these models are classified into discrete cell-based models and continuum models. At the tissue level, continuum models provide a very goodapproximation. These models incorporate a system of partial differential equations (PDEs), wherecell density, nutrients (i.e., oxygen and glucose), etc., are tracked. Modeling, mathematical analysisand numerical simulations were carried out in numerous papers, see [2, 5, 6, 7, 9, 11, 17, 21, 25,26, 30, 28] and the references therein. Lowengrub et al [19] provided a systematic review of tumormodel studies. However, in many of these models, the movement of the ECM within the tumor cellsis ignored. Therefore, in order to better describe and understand the whole process and relatedmechanism, we study a mathematical model for the influence of the extracellular matrix (ECM) on
Corresponding author: Rui Li.keywords: Free boundary, Shooting method, Singular integro-differential equations, Uniqueness, Simulation.2010 Mathematics Subject Classification: 35R35, 35M10, 35L45, 34B16. a r X i v : . [ m a t h . A P ] A ug RUI LI AND BEI HU tumor’s evolution in terms of system of partial differential equation. This model basically consistsof a system of parabolic equations and a hyperbolic equation for the density of the nutrient, forthe matrix degrading enzymes (MDE) and the ECM concentration. Moreover, our model is moreflexible, since it involves nonconstant coefficient µ ( E ) and allows the movement of the ECM fibres.All these considerations make our model into a more reasonable and realistic setting, but lead to amore challenging problem to analyze. 2. The model
In this section, we consider a PDE system to describe the evolution of the tumor.2.1.
Nutrients.
Let Ω( t ) denote the tumor domain at time t , and nutrient σ within the tumor is modeled by adiffusion equation c ∂σ∂t = ∆ σ + γ ( σ B − σ ) − λσ in Ω( t ) , (2.1)where c = T diffusion /T growth is the ratio of the nutrient diffusion time scale to the tumor growth(e.g. tumor doubling) time scale, γ ( σ B − σ ) denotes the nutrient supplied by the vasculature with γ being the transfer rate of nutrient in blood to tissue and σ B being the concentration of nutrientsin the vasculature, λσ describes the rate of consumption by the tumor. By appropriate change ofvariables, (2.1) is reduced to (see [14]) c ∂σ∂t = ∆ σ − λσ in Ω( t ) . (2.2)2.2. Extracellular Matrix.
The concentration of the ECM in the system is governed by contributions from three factors:haptotaxis, degrading, production. Here, there is a basic assumption that an equilibrium amountof nutrient ¯ σ is needed for tumor to sustain itself; beyond this ¯ σ the tumor grows, and below ¯ σ thetumor shrinks. Therefore a linear approximation for the proliferation S is given by S = µ ( σ − ¯ σ ) (¯ σ > , (2.3)where µσ represents the growth rate and µ ¯ σ represents the death rate from apoptosis. We shallemploy Darcy’s law (see [5, 6, 11]): −→ V = − ˜ µ ∇ p, (2.4)where −→ V represents the velocity of proliferating cells and p the pressure within the tumor resultingfrom this proliferation. It is well known that Darcy’s law describes the velocity of fluid in a porousmedium, with the coefficient ˜ µ depending on the density of the porous medium, representing amobility that reflects the combined effects of cell-cell and cell-matrix adhesion. In employingDarcy’s law, we have assumed ECM to be the porous medium; ˜ µ = ˜ µ ( E ) depends on the amountof ECM present in the tumor. By conservation of massdiv −→ V = S. (2.5)Substituting (2.3) into (2.5), we obtain div −→ V = µ ( σ − ¯ σ ) . (2.6)In papers ([5, 6, 11, 14]) both µ and ˜ µ are assumed to be constants; these are good approximationswhen ECM dose not vary much. Here, we shall incorporate a more reasonable assumption thatboth µ and ˜ µ also depend on ECM density E . It is clear that µ ( E ) and ˜ µ ( E ) are both monotonedecreasing functions bounded from above and below by positive constants. In order to do so, wealso need to incorporate the equation for E (see [7, 8, 16]): PARABOLIC-HYPERBOLIC SYSTEM MODELING 3 ∂E∂t + div( E · −→ V ) = − γmE + φ ( E ) , (2.7)where the term div( E · −→ V ) represents the movement of ECM owing to the cell proliferation −→ V ; theterm − γmE represents the degrading of ECM by MDE, here m represents the concentration ofMDE (see [25]); finally the term φ ( E ) is a positive term representing reorganization of ECM. Sincethe growth rate of ECM is smaller when the ECM is denser, φ ( E ) is a positive monotone decreasingfunction of E .2.3. Matrix degrading enzymes.
MDE is produced by the tumor to degrade ECM so that the cells can escape. The equation forMDE is given by (see [1, 7]) ∂m∂t = D m ∆ m + α − βm in Ω( t ) , (2.8)where ∆ m represents diffusion, D m is the constant diffusion coefficient and − βm represents naturaldecay. Here we assume α to be a constant production rate by the tumor.To summarise, the model studied in this paper is as follows: c ∂σ∂t = ∆ σ − λσ in Ω( t ) , (2.9a) ∂E∂t + div( E · −→ V ) = − γmE + φ ( E ) in Ω( t ) , (2.9b) ∂m∂t = D m (cid:52) m + α − βm in Ω( t ) , (2.9c) −→ V = − ˜ µ ( E ) ∇ p in Ω( t ) , (2.9d)div −→ V = µ ( E )( σ − ¯ σ ) in Ω( t ) . (2.9e)2.4. Boundary and initial conditions.
We impose boundary conditions σ = 1 on ∂ Ω( t ) , (2.10) ∂m∂n = 0 on ∂ Ω( t ) , (2.11) p = κ on ∂ Ω( t ) . (2.12)Equation (2.10) represents a condition that the tumor is immersed in an environment of constantnutrients; equation (2.11) represents no exchange of MDE on the tumor boundary; and equation(2.12) represents the cell-to-cell adhesiveness, where κ is the mean curvature. Finally, assumingthe velocity is continuous up to the boundary, then V n = −→ V · −→ n = − ˜ µ ∇ p · −→ n on ∂ Ω( t ) , (2.13)where V n represents the velocity of the boundary ∂ Ω( t ) in the normal direction.Initial conditions:Ω(0) = Ω , σ | t =0 = σ ( x ) , m | t =0 = m ( x ) , E | t =0 = E ( x ) . (2.14)In comparison with a system assuming ECM to be constants, our system is more reasonable andcomplex because we assume that ECM satisfies a hyperbolic equation coupled with nutrient σ andpressure p deriving from cells’ proliferation. In this paper we shall study the radially symmetriccase. While tumors in vivo are not spherical, tumors in vitro are typically of spherical shape[29]. The structure of this paper is as follows. In section 3, we proceed to derive estimates toestablish global existence and uniqueness and gave the lower bounds estimate of tumor radius RUI LI AND BEI HU R ( t ). In section 4, we prove that there exists a unique stationary solution by the shooting method.In section 5, the corresponding numerical simulation confirms the expected asymptotic stabilityin certain parameter range. In appendix, we proved that the well-posed theorem for the generalsingular integro-differential equation, which is a preliminary work for section 4.3. Time dependent solution
In this section we are concerned with the existence of radially symmetric solution.3.1.
Reformulation of the radially symmetric problem.
In order to prove the existence ofthe solution, for convenience, we do a reformulation for the radially symmetric problem.In radially symmetric case −→ V = x | x | u , we have by (2.6)div −→ V = u r + 2 r u = µ ( E )( σ − ¯ σ ) , div( E −→ V ) = E · div −→ V + ( ∇ E ) · −→ V = Eµ ( E )( σ − ¯ σ ) + uE r . Substituting this into (2.7), we obtain ∂E∂t + u ∂E∂r = Q ( σ, m, E ) , (3.1)where Q ( σ, m, E ) = − γmE + φ ( E ) − Eµ ( E )( σ − ¯ σ ) . (3.2)Furthermore, in the radially symmetric case, the mean curvature κ is a constant, and once σ and E are determined, one can uniquely solve p from the following linear elliptic equation: (cid:40) − div( µ ( E )) ∇ p ) = µ ( E )( σ − ¯ σ ) for | x | < R ( t ) ,p = κ = R ( t ) for | x | = R ( t ) . Therefore, we can drop the equation for p . In summary, the equations in the radially symmetriccase are: c ∂σ∂t = 1 r ∂∂r (cid:18) r ∂σ∂r (cid:19) − λσ, ≤ r < R ( t ) , t > , (3.3) ∂E∂t + u ∂E∂r = Q ( σ, m, E ) , ≤ r < R ( t ) , t > , (3.4) ∂m∂t = D m r ∂∂r (cid:18) r ∂m∂r (cid:19) + α − βm, ≤ r < R ( t ) , t > , (3.5) u r + 2 r u = µ ( E )( σ − ¯ σ ) , ≤ r < R ( t ) , t > . (3.6)The system is supplemented with boundary conditions σ ( R ( t ) , t ) = 1 , t > , (3.7) ∂m∂r ( R ( t ) , t ) = 0 , t > , (3.8)and free boundary condition (assuming continuity of velocity up to the boundary) dR ( t ) dt = u ( R ( t ) , t ) . (3.9)The system is also supplemented with initial conditions σ ( r,
0) = σ ( r ) , ≤ r ≤ R (0) , (3.10) m ( r,
0) = m ( r ) , ≤ r ≤ R (0) , (3.11) E ( r,
0) = E ( r ) , ≤ r ≤ R (0) . (3.12) PARABOLIC-HYPERBOLIC SYSTEM MODELING 5
By symmetry, ∂σ∂r (0 , t ) = 0 , ∂E∂r (0 , t ) = 0 , ∂m∂r (0 , t ) = 0 , u (0 , t ) = 0 . (3.13)We now use a change of variables that transform the free boundary into a fixed boundary: rR ( t ) → ˜ r, t → ˜ t, σ → ˜ σ, m → ˜ m, E → ˜ E, uR → ˜ u. For simplicity, we drop ” ∼ ” in our notation, and then σ ( r, t ) , m ( r, t ) , E ( r, t ) , u ( r, t ) satisfy cσ t = 1[ R ( t )] ( σ rr + 2 r σ r ) + R (cid:48) ( t ) R ( t ) rσ r − λσ, < r < , t > , (3.14a) m t = D m R ( t )] ( m rr + 2 r m r ) + R (cid:48) ( t ) R ( t ) rm r + α − βm, < r < , t > , (3.14b) E t + vE r = Q ( σ, m, E ) , < r < , t > , (3.14c) v ( r, t ) = u ( r, t ) − ru (1 , t ) , < r < , t > , (3.14d) u r + 2 r u = µ ( E )( σ − ¯ σ ) , < r < , t > , (3.14e) R t ( t ) = R ( t ) u (1 , t ) , t > , (3.14f) σ r (0 , t ) = 0 , σ (1 , t ) = 1 , t > , (3.14g) m r (0 , t ) = 0 , m r (1 , t ) = 0 , t > , (3.14h) u (0 , t ) = 0 , t > , (3.14i)with initial conditions R (0) = R , (3.15a) σ ( r,
0) = σ ( r ) , < r < , (3.15b) m ( r,
0) = m ( r ) , < r < , (3.15c) E ( r,
0) = E ( r ) , < r < , (3.15d)where σ , m ∈ C [0 ,
1] and E ∈ C [0 ,
1] satisfy R > , σ r (0) = 0 , σ (1) = 1 , m r (0) = 0 , m r (1) = 0 , E r (0) = 0 . (3.16)In particular, from (3.14i) and (3.14d), we see that v (0 , t ) = 0 , v (1 , t ) = 0 . This implies that no boundary conditions are needed for E at r = 0 , σ , m are radially symmetric functions and belong to W ,p ( B ) (for some fixed p > E ∈ C [0 , R >
0. By biological consideration, these initial functions are nonnegative and do notvanish completely.3.2.
Local existence and uniqueness.
We start with local existence and uniqueness by applyingthe contraction mapping principle.For
T >
0, we set Q T = (0 , × (0 , T ) , ¯ Q T = [0 , × [0 , T ] . Definition 3.1.
For any given
T > , we define a complete metric space ( X T , d ) as follows: X T is the subset of W ,p (0 , T ) × C , ( ¯ Q T ) consisting of a collection of pairs of functions R = R ( t ) , E = E ( r, t ) satisfying RUI LI AND BEI HU (i) R ∈ W ,p (0 , T ) , R (0) = R and (cid:107) R (cid:48) (cid:107) L p (0 ,T ) ≤ , R ≤ R ( t ) ≤ R , t ∈ [0 , T ] . (3.17) (ii) E ∈ C , ( ¯ Q T ) , E ( r,
0) = E ( r ) and (cid:107) E (cid:107) C , ( ¯ Q T ) ≤ M , (3.18) where M = (cid:107) E (cid:107) L ∞ + 1 .The metric d in X T is as follows d (( R , E ) , ( R , E ))= (cid:107) R − R (cid:107) C [0 ,T ] + (cid:107) R (cid:48) − R (cid:48) (cid:107) L p (0 ,T ) + (cid:107) E − E (cid:107) C , ( ¯ Q T ) . Given a pair (
R, E ) ∈ X T , we define σ = σ ( r, t ) as the solution of the initial boundary valueproblem (3.14a), (3.14g), (3.15b), and m = m ( r, t ) as the solution of the initial boundary valueproblem (3.14b), (3.14h), (3.15c). Let us define u = u ( r, t ) as the solution of (3.14e), (3.14i) and v = v ( r, t ) by (3.14d), that is u ( r, t ) = 1 r (cid:90) r g ( σ ( ρ, t ) , E ( ρ, t )) ρ dρ, (3.19)where g ( σ, E ) = µ ( E )( σ − ¯ σ ). We next define ( ˆ R, ˆ E ) = F ( R, E ) such thatˆ E t ( r, t ) + v ( r, t ) ˆ E r ( r, t ) = Q ( σ, m, ˆ E ) , < r < , t > , (3.20)ˆ E ( r,
0) = E ( r ) , < r < , (3.21)ˆ R (cid:48) = u (1 , t ) ˆ R ( t ) , t > , (3.22)ˆ R (0) = R . (3.23)We shall prove the existence of a local solution of (3.14)-(3.16) by using the contraction mappingtheorem for a map F : X T → X T .Clearly, the problem (3.22)-(3.23) can be solved explicitlyˆ R ( t ) = R exp (cid:18)(cid:90) t u (1 , τ ) dτ (cid:19) . (3.24)To uniquely solve (3.20)-(3.21), we introduce the characteristic curves ending at ( r, t ) dξ ( r, t ; s ) ds = v ( ξ ( r, t ; s ) , s ) ,ξ ( r, t ; t ) = r. (3.25)Since v ( r, t ) is continuous in ( r, t ) and Lipschitz in r , the characteristic curves ξ is well defined for0 ≤ s ≤ t. We then rewrite (3.20) in the form dds ˆ E ( ξ ( r, t ; s ) , s ) = Q ( σ ( ξ ( r, t ; s ) , s ) , m ( ξ ( r, t ; s ) , s ) , ˆ E ( ξ ( r, t ; s ) , s )) . Note that v satisfies (3.14d), we see that the characteristic curves do not leave and enter thespace interval (0 , Q ( r, t, E ) = Q ( σ ( r, t ) , m ( r, t ) , E ) andconsider d ˜ E ( r, t ; s ) ds = Q ( ξ ( r, t ; s ) , s, ˜ E ( r, t ; s )) , < s < t, ˜ E ( r, t ; 0) = E ( ξ ( r, t ; 0)) . (3.26)Clearly, (3.26) admits a unique (local) solution ˜ E . Thus, ˆ E ( r, t ) = ˜ E ( r, t ; t ) is the solution of(3.20)-(3.21). PARABOLIC-HYPERBOLIC SYSTEM MODELING 7
If we regard equation (3.14a) as a 1-dimensional parabolic equation with the spatial variable r ,then the coefficient of ∂σ/∂r has singularity at tumor center r = 0 due to∆ σ = ∂ σ∂r + 2 r ∂σ∂r . However, this singularity can be eliminated by employing the three-dimensional Cartesian coordi-nate.Due to the assumption imposed on R ( t ) in (3.17), we see that the coefficient R (cid:48) /R in equation(3.14a) and equation (3.14b) only belongs to L p . One can apply the classical parabolic theoryto obtain the strong solution σ and m exist and belong to W , ,p ( Q T ), see Theorem 9.1 and itscorollary of chapter IV in [20].In order to prove F maps X T into itself for some small T , it suffices to estimate the norms of( ˆ R, ˆ E ) as well as σ, m, u, v .For notational convenience, in the sequel we shall denote by C any one of several constants whichdepend on R , M , but does not depend on T ∈ (0 , Lemma 3.2. If R ( t ) satisfies (3.17) , then the strong solutions σ and m admit the following uniformbounds (cid:107) σ (cid:107) W , ,p ( B × (0 ,T )) ≤ C, (cid:107) σ (cid:107) C α, (1+ α ) / ( ¯ B × [0 ,T ]) ≤ C, (3.27) (cid:107) m (cid:107) W , ,p ( B × (0 ,T )) ≤ C, (cid:107) m (cid:107) C α, (1+ α ) / ( ¯ B × [0 ,T ]) ≤ C, (3.28) where α = 1 − /p and B is the unit ball in R . Proof : The bounds for σ and m are similar, we focus the uniform estimates for σ only. Note thatall functions are defined in the time interval [0 , T ], we can extend R ( t ), so that it is defined in afixed interval [0,1]. More precisely, ˜ R ∈ W ,p (0 ,
1) is defined as follows˜ R ( t ) = (cid:40) R ( t ) , for t ∈ [0 , T ] ,R ( T ) , for t ∈ ( T, . It is clear that ˜ R (cid:48) ( t ) ≡ T < t < . Clearly, (cid:107) ˜ R (cid:107) W ,p (0 , ≤ (cid:107) ˜ R (cid:107) C [0 , + (cid:107) ˜ R (cid:48) (cid:107) L p (0 , ≤ R + 1 . Since the embedding W ,p (0 , (cid:44) → C α [0 ,
1] is continuous, we conclude that 1 / ˜ R is continuous and (cid:107) ˜ R (cid:107) C α [0 , ≤ C (cid:107) ˜ R (cid:107) W ,p (0 , ≤ C. We now define ˜ σ to be the solution of (3.14a), (3.14g), (3.15b) (with R replaced by ˜ R ) in the timeinterval [0 , L p theory [20], we see that ˜ σ ∈ W , ,p ( B × (0 , T )) exists and (cid:107) ˜ σ (cid:107) W , ,p ( B × (0 , ≤ C. (3.29)Recalling that the embedding W , ,p ( B × (0 , (cid:44) → C α, (1+ α ) / ( ¯ B × [0 , (cid:107) ˜ σ (cid:107) C α, (1+ α ) / ( ¯ B × [0 , ≤ C (cid:107) ˜ σ (cid:107) W , ,p ( B × (0 , ≤ C. (3.30)By uniqueness of parabolic equation, we see ˜ σ and σ are the same in the time interval [0 , T ). Theconclusion immediately follows from (3.29) and (3.30). RUI LI AND BEI HU
From (3.19), we deduce that (cid:107) u (cid:107) C , ( ¯ Q T ) ≤ C, (cid:107) v (cid:107) C , ( ¯ Q T ) ≤ C. (3.31)From (3.24) and (3.22), we obtain that (cid:107) ˆ R (cid:107) C [0 ,T ] ≤ C. (3.32)In particular, for sufficiently small T > t ∈ [0 ,T ] | ˆ R ( t ) − R | ≤ CT ≤ R / , | ˆ R (cid:48) ( t ) | L p (0 ,T ) ≤ CT /p ≤ . (3.33)Therefore, ˆ R = ˆ R ( t ) satisfies (3.17) in the definition of X T provided T is sufficiently small. Lemma 3.3.
For any sufficiently small
T > , if ( R, E ) ∈ X T , then the unique solution ˆ E of (3.20) - (3.21) is well-defined for t ∈ [0 , T ] and satisfies (cid:107) ˆ E (cid:107) C , ( ¯ Q T ) ≤ C. (3.34) Proof : Note that 0 ≤ ˜ E ( r, t ; 0) ≤ (cid:107) E (cid:107) ≤ M and Q ( r, t, ≥ , Q ( r, t, E ) ≤ Ψ( E )for r ∈ [0 , t ∈ (0 , T ], E ≥
0, where Ψ( E ) is some positive smooth function defined in [0 , ∞ ). Let y ( s ) denote the unique solution of the initial problem dyds = Ψ( y ) , y (0) = M which exists at least in a finite interval [0 , h ] for some h >
0. Then by the comparison principle forODE, we deduce that the solution ˜ E of (3.26) satisfies0 ≤ ˜ E ( r, t ; s ) ≤ y ( s ) ,s ∈ [0 , h ], where ˜ E is a solution of (3.26). This shows that ˜ E exists and is bounded in [0 , × [0 , h ] × [0 , h ] provided T > h . In particular, ˆ E exist and is bounded in ¯ Q T for T ≤ h .From (3.31) and (3.14d), we see that (cid:107) v (cid:107) C , ≤ C . It follows that ξ , defined by (3.25), belongsto C , , and satisfies (cid:107) ξ (cid:107) L ∞ + (cid:107) ξ r (cid:107) L ∞ + (cid:107) ξ t (cid:107) L ∞ + (cid:107) ξ s (cid:107) L ∞ ≤ C. Recalling that Q ( r, t, E ) is C α in r , C (1+ α ) / in t and smooth in ˜ E , we deduce from (3.26) that | ˜ E r | < C and hence | ˆ E r | < C in ¯ Q T , here we write ˆ E r ( r, t ) = ˜ E r ( r, t ; t ). Finally, we derive from equation (3.14c) that ˆ E t ( r, t ) = − v ( r, t ) ˆ E r ( r, t ) + Q ( σ, m, ˆ E ) is also bounded in ¯ Q T . This completes the proof. Corollary 3.4.
For sufficiently small
T > , if ( R, E ) ∈ X T , then ˆ E satisfies (3.18) . Proof : From Lemma 3.3, we have | ˆ E ( r, t ) | ≤ (cid:107) E (cid:107) L ∞ + CT.
This shows that ˆ E = ˆ E ( r, t ) satisfies (3.18) provided T is sufficiently small.Combing (3.33) and Corollary 3.4, we have established the following Proposition 3.5.
For sufficiently small
T > , F is well-defined and maps X T into itself. PARABOLIC-HYPERBOLIC SYSTEM MODELING 9
We now establish that F is a contraction mapping for sufficiently small T . Let ( R i , E i ) ∈ X T for i = 1 ,
2. Let σ i , m i , u i , v i , ( ˆ R i , ˆ E i ) = F ( R i , E i ) be the solution corresponding to ( R i , E i ). Recall d = d (( R , E ) , ( R , E ))= (cid:107) R − R (cid:107) C [0 ,T ] + (cid:107) R (cid:48) − R (cid:48) (cid:107) L p (0 ,T ) + (cid:107) E − E (cid:107) C , ( ¯ Q T ) . (3.35)From equation (3.14a), we see σ ∗ = σ − σ satisfies c ∂σ ∗ ∂t = R ( t ) ( ∂ σ ∗ ∂r + r ∂σ ∗ ∂r ) + R (cid:48) ( t ) R ( t ) r ∂σ ∗ ∂r − λσ ∗ + h ( r, t ) , ( σ ∗ ) r (0 , t ) = 0 , σ ∗ (1 , t ) = 0 ,σ ∗ ( r,
0) = 0 , where h ( r, t ) = ( 1 R − R )( ∂ σ ∂r + 2 r ∂σ ∂r ) + ( R (cid:48) R − R (cid:48) R ) r ∂σ ∂r . From the estimates of σ (Lemma 3.2), we derive (cid:107) h (cid:107) L p ( B × (0 ,T )) ≤ C ( (cid:107) R − R (cid:107) C [0 ,T ] + (cid:107) R (cid:48) − R (cid:48) (cid:107) L p (0 ,T ) ) . Using the same method as that in Lemma 3.2, we deduce (cid:107) σ − σ (cid:107) W , ,p ( B × (0 ,T )) ≤ Cd, (cid:107) σ − σ (cid:107) C α, (1+ α ) / ( ¯ B × [0 ,T ]) ≤ Cd. (3.36)Similarly, (cid:107) m − m (cid:107) W , ,p ( B × (0 ,T ) ≤ Cd, (cid:107) m − m (cid:107) C α, (1+ α ) / ( ¯ B × [0 ,T ]) ≤ Cd. (3.37)By the definition of u i and v i we have (cid:107) u − u (cid:107) C , ( ¯ Q T ) ≤ C ( (cid:107) σ − σ (cid:107) ∞ + (cid:107) E − E (cid:107) ∞ ) ≤ Cd, (cid:107) v − v (cid:107) C , ( ¯ Q T ) ≤ Cd. (3.38)Set R ∗ = ˆ R − ˆ R . Then, by direct calculations we see that R ∗ satisfies (cid:40) dR ∗ ( t ) dt = ¯ g ( t ) R ∗ ( t ) + ¯ g ( t ) ,R ∗ (0) = 0 , (3.39)or R ∗ ( t ) = (cid:90) t ¯ g ( τ ) e (cid:82) tτ g (˜ τ ) d ˜ τ dτ, (3.40)where ¯ g ( t ) = u (1 , t ), ¯ g ( t ) = R ( t )( u (1 , t ) − u (1 , t )) satisfy | ¯ g ( t ) | ≤ C, | ¯ g ( t ) | ≤ Cd, t ∈ [0 , T ] (3.41)by using (3.31), (3.32) and (3.38). One can easily derive from (3.39)-(3.41) that (cid:107) ˆ R − ˆ R (cid:107) C [0 ,T ] ≤ CT d, (cid:107) ˆ R (cid:48) − ˆ R (cid:48) (cid:107) C [0 ,T ] ≤ Cd, (cid:107) ˆ R (cid:48) − ˆ R (cid:48) (cid:107) L p (0 ,T ) ≤ CT /p d. (3.42) Lemma 3.6.
There holds (cid:107) ˆ E − ˆ E (cid:107) C ( ¯ Q T ) ≤ CT d (3.43) for sufficiently small T . Proof : Note that ˆ E ∗ = ˆ E − ˆ E satisfies ∂ ˆ E ∗ ∂t + ∂ ˆ E ∗ ∂r v − A ( r, t ) ˆ E ∗ = h ( r, t ) , ˆ E ∗ ( r,
0) = 0 , (3.44)where A ( r, t ) = (cid:90) Q ( σ ( r, t ) , m ( r, t ) , θ ˆ E ( r, t ) + (1 − θ ) ˆ E ( r, t )) dθ,h ( r, t ) = − ( v − v ) ∂ ˆ E ∂r + [ Q ( σ , m , ˆ E ) − Q ( σ , m , ˆ E )]+ [ Q ( σ , m , ˆ E ) − Q ( σ , m , ˆ E )] . Using the estimates (3.27),(3.36) for σ i , the estimates (3.28),(3.37) for m i , the estimates (3.31),(3.38)for u i , v i and the estimates (3.34) for ˆ E i , we get | A ( r, t ) | ≤ C, | h ( r, t ) | ≤ Cd, ( r, t ) ∈ ¯ Q T . Hence, integrating (3.44) along its characteristics as before, we find that (cid:107) ˆ E ∗ (cid:107) C ( ¯ Q T ) ≤ CT d.
This finishes the proof.Now we show the local existence and uniqueness of solution to (3.14)-(3.15).
Theorem 3.7.
Assume that the initial data satisfy (3.16) . Then there exists a
T > which onlydepends on (cid:107) σ (cid:107) W ,p ( B ) , (cid:107) m (cid:107) W ,p ( B ) , R , and (cid:107) E (cid:107) C [0 , , such that problem (3.14) - (3.15) admitsa unique solution σ, m ∈ W , ,p ( B × (0 , T )) , u, v ∈ C , ( ¯ Q T ) , E ∈ C , ( ¯ Q T ) , R ∈ C [0 , T ] . Proof : From Proposition 3.5, (3.42) and Lemma 3.6, we deduce that for sufficiently small
T > F is a contraction mapping from X T into itself. Therefore, there is a unique fixed point in X T for F , and thus (3.14)-(3.15) admits a unique solution ( σ, m, E, u, v, R ) in the time interval [0 , T ]. Remark 3.8.
This theorem shows the local existence and uniqueness of radially symmetric solutionof our equation. If the E satisfies ( E ) r (0) = 0 , then one can prove that E r (0 , t ) = 0 for all t .If initial data σ , m , E are radially symmetric function, and σ , m ∈ C ,α ( ¯ B ) (for some α ∈ (0 , ), E ∈ C ( ¯ B ) , R > , then the standard PDE theory tells us the solution obtained inTheorem 3.7 are also radially symmetric and are more regular, i.e. σ, m ∈ C α, (2+ α ) / ( ¯ B × [0 , T ]) , u, v ∈ C , ( ¯ Q T ) , E ∈ C ( ¯ B × [0 , T ]) , R ∈ C [0 , T ] . Global existence.Theorem 3.9.
Assume that the initial data satisfy (3.16) . Then problem (3.14) - (3.15) admits aunique solution ( σ, m, u, v, E, R ) globally in time. Proof : Suppose to the contrary that [0 , ˜ T ) is the maximum time interval ( ˜ T < ∞ ) for the existenceof the solution. We proceed to derive necessary estimates for the global existence.Employing the maximum principle for parabolic equations, we deduce that0 < σ ( r, t ) , m ( r, t ) ≤ max {(cid:107) σ (cid:107) L ∞ , αβ , (cid:107) m (cid:107) L ∞ } (3.45)for all ( r, t ) ∈ [0 , × [0 , ˜ T ). We shall denote by C various constant which is independent of ˜ T , andby C ( ˜ T ) various constants which depends on ˜ T ∈ (0 , ∞ ). Since µ ( E ) is decreasing and bounded,we obtain from (3.19) and (3.45) that | u ( r, t ) | ≤ C, | u r ( r, t ) | ≤ C for all ( r, t ) ∈ [0 , × [0 , ˜ T ) . (3.46) PARABOLIC-HYPERBOLIC SYSTEM MODELING 11
From equation (3.14f) and (3.46), we derive that R e − C ˜ T ≤ R ( t ) ≤ R e C ˜ T , | R (cid:48) ( t ) | ≤ CR e C ˜ T for all t ∈ [0 , ˜ T ) . (3.47)By the L p theory (see [20] or [18]) and Sobolev inequalities, we get (cid:107) σ ( r, t ) , m ( r, t ) (cid:107) W , ,p ( B × [0 , ˜ T )) ≤ C ( ˜ T ) , (cid:107) σ ( r, t ) , m ( r, t ) (cid:107) C α, α ( ¯ B × [0 , ˜ T ]) ≤ C ( ˜ T ) . (3.48)From the bounds (3.45) for σ and m and the assumption on Q , we obtain that Q ( r, t, E ) = Q ( σ ( r, t ) , m ( r, t ) , E ) satisfies |Q ( r, t, E ) | ≤ C ∗ ( E + 1) , r ∈ [0 , , t ∈ [0 , ˜ T ) , E ≥ C ∗ which is independent of ˜ T ∈ (0 , ∞ ). Note that E ( r, t ) = ˜ E ( r, t ; t ) and ˜ E ( r, t ; s )satisfies d ˜ E ( r, t ; s ) ds = Q ( ξ ( r, t ; s ) , s, ˜ E ( r, t ; s )) , < s < t, ˜ E ( r, t ; 0) = E ( ξ ( r, t ; 0)) . (3.49)We conclude that | ˜ E ( r, t ; s ) | ≤ ( (cid:107) E (cid:107) C [0 , + 1) e C ∗ ˜ T − , r ∈ [0 , , ≤ s ≤ t < T. By using the bounds (3.46) for u and then for v , we see that the character curves ξ ( r, t ; s ) and itsderivatives ξ r ( r, t ; s ) are bounded by some constant C ( ˜ T ). Combining this with (3.50), one obtain | ˜ E r ( r, t ; s ) | ≤ C ( ˜ T ) , r ∈ [0 , , ≤ s ≤ t < T. This gives the bounds for E r and hence for E t by equation (3.14c). Therefore, | E ( r, t ) | + | E r ( r, t ) | + | E t ( r, t ) | ≤ C ( ˜ T ) for all ( r, t ) ∈ [0 , × [0 , ˜ T ) . (3.50)Taking ˜ T − (cid:15) (where 0 < (cid:15) < ˜ T is arbitrary) as a new initial time, then we can extend thesolution to Q ( ˜ T − (cid:15) )+ δ for some small δ > δ depends only on an upper bound of the data at time ˜ T − (cid:15) and the the lower bounds of R ( ˜ T − (cid:15) ). By a priori estimate (3.46)-(3.50), we find that δ dependsonly on ˜ T (but δ is independent of (cid:15) ), i.e., δ = δ ( ˜ T ). If we take (cid:15) < δ ( ˜ T ), then we get( ˜ T − (cid:15) ) + δ > ˜ T , which contradicts the assumption that [0 , ˜ T ) is the maximum time interval for the existence of thesolution. Therefore, the maximum time interval for the existence of the solution is [0 , ∞ ).3.4. The lower bound of R ( t ) . In this subsection, we consider the case of µ ( E ) is near 0 andstudy the lower bound of R ( t ).From (3.3), we get u ( t ) = 1 R ( t ) (cid:90) R ( t )0 µ ( E )( σ − ¯ σ ) ρ dρ. Combining with equation (3.9), we have R R (cid:48) ( t ) = (cid:90) R ( t )0 µ ( E )( σ − ¯ σ ) ρ dρ. (3.51)From the maximum principle for parabolic equation, we see that0 < σ ( r, t ) < ˆ σ for 0 ≤ r ≤ R ( t ) , t > , where ˆ σ := max { , max σ ( r ) } . Note that µ ( E ) is bounded, η < µ ( E ( r, t )) < η , for 0 ≤ r ≤ R ( t ) , t > for some 0 < η < η . Theorem 3.10.
There exists a positive constant δ , such that R ( t ) ≥ δ for all t > . Proof : We will prove this theorem in two steps:
Step 1.
We claim that lim sup t →∞ R ( t ) > . We shall argue by contradiction. If the conclusionis false, then lim t →∞ R ( t ) = 0 . (3.53)As in the reference [15], we let v ( r, t ) = ¯ σ R ( t ) r sinh( M r )sinh(
M R ( t )) for t ≥ t , (3.54)where M = λ + 2 + N and N is a positive number. As r → + , we have that r sinh r = 1 − r + 7360 r + O ( r ) , sinh rr ddr (cid:18) r sinh r (cid:19) = − r + 145 r + O ( r ) . Note µ ( E ) , σ are bounded, we assume that | µ ( E )( σ − ¯ σ ) | < η, r < R ( t ) , t > . Take a small δ satisfying 0 < δ < (cid:114) M η , then for R ∈ (0 , δ ) − M R < (cid:18) sinh rr ddr (cid:18) r sinh r (cid:19)(cid:19)(cid:12)(cid:12)(cid:12)(cid:12) r = MR < , > M R sinh( M R ) > max (cid:26) ηM δ , σ (cid:27) . From (3.51), | R (cid:48) ( t ) | < ηR ( t )3 for all t > . Thus, (cid:12)(cid:12)(cid:12)(cid:12) dvdt (cid:12)(cid:12)(cid:12)(cid:12) = (cid:12)(cid:12)(cid:12)(cid:12) v ( r, t ) sinh rr ddr (cid:18) r sinh r (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) r = MR | M R (cid:48) ( t ) |≤ M R ( t ) | R (cid:48) ( t ) |≤ ηM R ( t )) . Then, if R ( t ) < δ , we have v t − ∆ v + λv = v t − v − N v< v t − v ≤ ηM R − ηM δ < . From (3.53), there is a large time T such that R ( t ) < δ , t ≥ T and v t − ∆ v + λv ≤ , if r < R ( t ) , t ≥ T . (3.55) PARABOLIC-HYPERBOLIC SYSTEM MODELING 13
Consider the function w = σ − v + z, where z = e − λ ( t − T ) . It satisfies w t − ∆ w + λw ≥ , if r < R ( t ) , t ≥ T , and it is positive on r = R ( t ) , t > T and on { t = T , r < R ( T ) } . By the maximum principle, w > t > t , i.e., σ ( r, t ) ≥ v ( r, t ) − z ( t ) . This inequality can be used to estimate R (cid:48) from below: R ( t ) R (cid:48) ( t ) = (cid:90) R ( t )0 µ ( E )( σ ( r, t ) − ¯ σ ) ρ dρ ≥ (cid:90) R ( t )0 µ ( E )( v − ¯ σ ) ρ dρ − (cid:90) R ( t )0 µ ( E ) z ( t ) ρ dρ ≥ η − ¯ σ ) R ( t ) − η e − λ ( t − T ) R ( t ) > . if t is sufficiently large by (3.52). It follows that for some large time T ( T > T ) R (cid:48) ( t ) > t > T , i.e. R ( t ) is monotone increasing. This contradicts to (3.53). Thus we finish the proof of step 1. Step 2.
We claim that lim inf t −→∞ R ( t ) ≥ δ , where δ = θδ and δ = θ δ . To prove the aboveresult, suppose for contradiction that lim inf t −→∞ R ( t ) < δ . Then by the step 1, there exists a t = t such that R ( t ) = δ . We shall prove that R ( t ) ≥ δ for all t > t , (3.56)and this establishes the theorem. Suppose that (3.56) is false, then there exists t < t < t suchthat R ( t ) = δ , R ( t ) = δ and δ < R ( t ) < δ for t < t < t ,R (cid:48) ( t ) ≤ . (3.57)By (3.14f) and R (cid:48) /R ≥ − η , we have that t − t ≥ η log δ δ = 1 η log 1 θ = γ,t − t ≥ η log δ δ = 1 η log 1 θ = γ. The domain D = { ( r, t ) : r < ρ = δ e ηγ , t − γ < t < t } contains the sub-domain D = { ( r, t ) : r < R ( t ) , t − γ < t < t } . We introduce the solution W to W t = ∆ W − λW, in D ,W = 1 , r = ρ , t ∈ ( t − γ, t ) ,W = 0 , r < ρ , t = t − γ. (3.58)By comparison with Green’s function for a rectangular domain constructed by a series of reflections,we obtain W ( r, t ) ≥ (cid:15) > , r < ρ , then by maximum principle, we have σ ≥ w in D , thus σ ( r, t ) ≥ (cid:15) , r < ρ . Next we introduce the domain D = { ( r, t ) : r < R ( t ) , t < t < t } (3.59)and a comparison function in D : v ( r, t ) = ˆ σ ( t ) R ( t ) r sinh( M r )sinh(
M R ( t )) , where ˆ σ ( t ) = (cid:15) e N ( t − t ) , t < t < t , (3.60) N = 1 t − t log 1 (cid:15) , so that ˆ σ ( t ) = (cid:15) , ˆ σ ( t ) = 1 . As in Step 1 we compute v t − ∆ v + λv = ˆ σ (cid:48) ( t )ˆ σ ( t ) v + ˆ σ ( t )( ∂ t − ∆ + λ ) v ˆ σ ( t )= N v + ˆ σ ∂∂t (cid:18) v ˆ σ ( t ) (cid:19) − v − N v = ˆ σ ∂∂t ( v ˆ σ ( t ) ) − v ≤ ˆ σ ( t ) (cid:20) ηM R − o ( R )) (cid:21) < . Thus v is a subsolution. In view of (3.59) and (3.60), we see that σ ≥ v on both r = R ( t ) and t = t . Here the maximum principle implies σ > v in D and σ ( r, t ) ≥ R ( t ) r sinh( M r )sinh(
M R ( t ))and σ ( r, t ) ≥ M δ sinh( M δ ) > σ . Using this in (3.51) we deduce that R (cid:48) ( t ) > , a contradiction to (3.57).4. The existence of radially symmetric stationary solution
In this section, we derive the existence and uniqueness of the radially symmetric stationarysolution. The major challenge for establishing existence and uniqueness stems from the singularityof our integro-differential equation. These types of equations are not covered by the standardtheory. Another challenge is to establish continuity of the velocity field near both ends r = 0 and r = R s . In addition, for such a highly nonlinear system, the uniqueness is by no means trivial. Asa matter of fact, uniqueness may not be valid for certain system (e.g., stationary problem for theprotocell [13]). We have explored the special structure of our problem which enables us to overcomethe difficulties and established uniqueness.In order to establish the existence, another important work in our paper is to construct thewell-posed of the general singular equations. Since the proof is lengthy and complex, we put it inthe appendix. PARABOLIC-HYPERBOLIC SYSTEM MODELING 15
We consider the general singular integro-equation (cid:40) dxdt = f ( x,t ) (cid:82) t g ( x ( s ) ,s ) ds ,x (0) = x , (4.1)where f , g are two functions defined in a domain G ⊂ R × R . We assume that f , g satisfy(F1) ( x , ∈ G and f ( x ,
0) = 0;(F2) f and its derivatives f x , f t are continuous in G ;(F3) g ∈ C ( G ) and g is local Lipschitz continuous with respect to x ;(F4) g ( x , (cid:54) = 0 and θ = f x ( x , /g ( x , ∈ ( −∞ , Theorem 4.1.
Assume f and g satisfy conditions (F1) though (F4). Then (4.1) admits a unique C solution in [ −T , T ] for a small T . Moreover, x (cid:48) (0) = f t ( x , g ( x , − f x ( x , . Remark 4.2.
The above theorem guarantees uniqueness in the class C [0 , T ] . In the case θ < ,the continuous solution x ∈ C [0 , T ] ∩ C (0 , T ] of (4.1) exists and is unique, and it is also in C class. However, it is crucial to notice that in the case θ ∈ (0 , , the continuous solution of (4.1) may not be unique, but the C solution is unique. Remark 4.3.
In fact, Theorem 4.1 tells us the solution exists in a small interval. If one denotesby y = y ( t ) the denominator of (4.1) , then x = x ( t ) , y = y ( t ) satisfy a integro-differential equation,and the solution can be expanded to a maximal existence interval. We next consider an equation of the type (4.1) involving a parameter µ : (cid:40) dxdt = f ( x,t,µ ) (cid:82) t g ( x ( s ) ,s,µ ) ds ,x (0) = ϕ ( µ ) , (4.2)where f and g are defined in an open set G × U of ( R × R ) × R m with ( x , ∈ G , µ ∈ U . Weassume that f and g satisfy:(F5) ϕ ( µ ) is a continuous function in a neighborhood of µ , and ϕ ( µ ) = x , f (0 , ϕ ( µ ) , µ ) ≡ f, g ∈ C ( G × U ) are local Lipschitz continuous with respect to x in G × U .(F7) The derivatives f x , f t exist and are continuous in a neighborhood of ( x , , µ ).(F8) g (cid:54) = 0 and f x /g < x, t, µ ) = ( x , , µ ). Theorem 4.4.
Let f, g be the functions satisfying conditions (F5)-(F8). Assume that the uniquesolution x = x ( t, µ ) of (4.2) exists in a bounded closed interval [ a, b ] with a < < b . Then thereexists a δ > such that(1) For every µ ∈ B µ ,δ , the solution x = x ( t, µ ) of (4.2) exists in [ a, b ] .(2) The function ( t, µ ) (cid:55)→ x ( t, µ ) is continuous in [ a, b ] × B µ ,δ . Next we consider steady state equations in the radially symmetric case are: r ∂∂r (cid:18) r ∂σ∂r (cid:19) − λσ = 0 , < r < R, (4.3a) D m r ∂∂r (cid:18) r ∂m∂r (cid:19) + α − βm = 0 , < r < R, (4.3b) uE (cid:48) = Q ( σ, m, E ) := − γmE + φ ( E ) − Eµ ( E )( σ − ¯ σ ) , < r < R, (4.3c) u (cid:48) + 2 r u = µ ( E )( σ − ¯ σ ) , < r < R, (4.3d) σ (cid:48) (0) = 0 , σ ( R ) = 1 , (4.3e) m (cid:48) (0) = 0 , m (cid:48) ( R ) = 0 , (4.3f) u (0) = 0 , u ( R ) = 0 . (4.3g)It is clear from (4.3b) and (4.3f) that m = αβ . The stationary solution for σ . We impose the following structural conditions: for some
N > ∂Q ( σ, m, E ) ∂E + µ ( E )¯ σ < , < σ ≤ , m = αβ , < E < ∞ ,Q ( σ, m, N ) < , < σ ≤ , m = αβ . (4.4)Condition (4.4) also implies that Q ( σ, m, E ) is decreasing in E >
0, that is, ∂Q ( σ, m, E ) ∂E < . Lemma 4.5.
Assume that there is a solution ( σ, m, E, u, R ) to (4.3) . Then σ ( r ) = R sinh( √ λR ) sinh( √ λr ) r , m = αβ . (4.5) u satisfies u (0) = 0 , u ( R ) = 0 , u ( r ) < for r ∈ (0 , R ) and u ( r ) = 1 r (cid:90) r µ ( E ( ρ ))( σ ( ρ ) − ¯ σ ) ρ dρ. (4.6) E ∈ C [0 , R ] satisfies the following singular equation E (cid:48) ( r ) = Q ( σ ( r ) , m, E ( r )) r (cid:82) r µ ( E ( ρ ))( σ ( ρ ) − ¯ σ ) ρ dρ , (4.7) and E (0) , E ( R ) satisfy Q ( σ ( r ) , m, E ( r )) = 0 at r = 0 , R. (4.8) Finally,
R > satisfies σ (0) = √ λR sinh( √ λR ) < ¯ σ. (4.9) Proof : We immediately obtain σ , given explicity by (4.5). Note that from (4.3c) and (4.3g), weget the C solution E satisfies Q ( σ ( r ) , m, E ( r )) = 0 at r = 0 , R. From (4.3d) with (4.3g) and the fact that ¯ σ ∈ (0 , u is represented as u ( r ) = 1 r (cid:90) r µ ( E ( ρ ))( σ ( ρ ) − ¯ σ ) ρ dρ or u ( r ) = 1 r (cid:90) rR µ ( E ( ρ ))( σ ( ρ ) − ¯ σ ) ρ dρ. PARABOLIC-HYPERBOLIC SYSTEM MODELING 17 (4.7) then is a re-statement of (4.3c).Note that u (0) = u ( R ) = 0, µ ( E ) is positive, and σ ( r ) − ˜ σ is strictly increasing, we derive that σ − ˜ σ admits exactly one interior root r ∈ (0 , R ) and R satisfies (4.9). Therefore, u is negative in(0 , R ).Since σ ( r ) also depends on R , we write it as σ ( r, R ). From the assumption (4.4) and the fact Q ( σ, m,
0) = φ (0) >
0, there exist a unique h = h ( r, R ) > Q ( σ ( r, R ) , m, h ) = 0 . (4.10)Moreover, the implicit function theorem implies that h ( r, R ) is a smooth function in two variable( r, R ) and that ∂h∂r = − Q σ ( σ ( r, R ) , m, h ( r, R )) Q E ( σ ( r, R ) , m, h ( r, R )) ∂σ ( r, R ) ∂r . (4.11)From (4.4) and Q σ = − µ ( E ) E <
0, we conclude that h is a decreasing function in r , h r ( r, R ) < , r > . We also have h is uniformly bounded by positive constants from below and above, N < h ( r, R ) < N, r ∈ [0 , R ] , (4.12)where N = h ( R, R ) > R (since σ ( R, R ) ≡ Q ( σ ( r, R ) , m, ξ ) = < ξ > h ( r, R ) , = 0 for ξ = h ( r, R ) ,> ξ < h ( r, R ) . (4.13)In order to find a solution ( σ, m, E, u, R ) of (4.3), it is equivalent to obtain a solution ( E, u, R )of the following equations uE (cid:48) = Q ( σ ( r, R ) , m, E ) , (4.14a) u (cid:48) + 2 r u = µ ( E )( σ ( r, R ) − ¯ σ ) , (4.14b) u ( R ) = 0 , E ( R ) = h ( R, R ) , (4.14c) u (0) = 0 , E (0) = h (0 , R ) , (4.14d)where σ ( r, R ) = R sinh( √ λR ) sinh( √ λr ) r . (4.15)Note that problem (4.14) is an ordinary differential system, and there are singularities at r = 0 , R .Since problem (4.14) is a boundary value problem, we will solve (4.14) by the shooting method.For any fixed R >
0, we solve the initial problem (
E, u ) of (4.14a)-(4.14c) near r = R , then find asuitable R such that ( E, u ) exists in [0 , R ] and satisfies condition (4.14d) at the other end.4.2.
The existence for stationary solution E .Solutions near r = R , i.e., < R − r (cid:28) . If we integrate (4.14b) with u ( R ) = 0, we get u ( r ) = 1 r (cid:90) rR µ ( E )( σ − ¯ σ ) ρ dρ. Therefore, E satisfies the following initial value problem: (cid:40) dE ( r ) dr = Q ( σ ( r ) ,m,E ( r )) r (cid:82) rR µ ( E ( ρ ))( σ ( ρ ) − ¯ σ ) ρ dρ ,E ( R ) = h ( R ) . (4.16) Lemma 4.6.
The equation (4.16) admits a unique solution. Moreover, dEdr ( R ) = E ( R ) µ ( E ( R )) σ r ( R ) ∂Q∂E (1 , m, E ( R )) − µ ( E ( R ))(1 − ¯ σ ) . (4.17) Proof : Note that E ( R ) = h ( R ) > Q E <
0, we derive˜ θ (cid:44) γm − φ E ( R ) − ∂ ( Eµ ( E )) ∂E (¯ σ − σ ) µ ( E ( R ))(¯ σ −
1) = − ∂Q∂E | r = R µ ( E ( R ))(¯ σ − < . Therefore, the conditions in Theorem 4.1 are satisfied. We can, after a change of variables r → R − r, apply Theorem 4.1 to show the existence and uniqueness of the solution E ∈ C [ R − δ, R ] to (4.3),for a small δ > . Lemma 4.6 tells us the solution E of (4.16) exists in the interval [ R − δ, R ], it also exists in aright neighborhood of r = R . The equations (4.16) or (4.14a)-(4.14b) do not have a singularityat r = R − δ . The standard ODE theory implies that the solution pair ( E, u ) can be extendedto a maximal interval ( τ, R ], so long as the denominator u is negative and E is positive on ( τ, R ).Noting that there is a singularity terms 2 u/r in equation (4.14b), we see that τ ≥
0. It is clear that E ( r ) > , u ( r ) < , r ∈ ( τ, R ) . (4.18)Since R is also unknown, we usually denote the solution by E ( r, R ), u ( r, R ) and denote by τ = τ ( R )for the maximum existence interval ( τ ( R ) , R ]. Remark 4.7.
The initial solution ( E, u ) of (4.14a) , (4.14b) and (4.14c) is also depending contin-uously on R ; see Theorem 4.4. More precisely, if for any fixed ¯ R > , the solution E ( r, ¯ R ) , u ( r, ¯ R ) of (4.14a) , (4.14b) and (4.14c) are well-solved in a compact interval [ˆ r, ˇ r ] with ¯ R ∈ (ˆ r, ˇ r ) , thenthere exists δ > such that(i) E ( r, R ) , u ( r, R ) of (4.14a) , (4.14b) and (4.14c) are well-solved in the interval [ˆ r, ˇ r ] for every R ∈ [ ¯ R − δ, ¯ R + δ ] .(ii) E ( r, R ) , u ( r, R ) are continuous functions in the variable ( r, R ) ∈ [ˆ r, ˇ r ] × [ ¯ R − δ, ¯ R + δ ] . Solution near < r (cid:28) . We consider dEdr = Q ( σ ( r ) , m, E ) r (cid:82) r µ ( E )( σ − ¯ σ ) ρ dρ . The singularity at r = 0 is not exactly of the format in Theorem 4.1, but the quantity θ defined in(F4) proceeding Theorem 4.1 is very close to˜ θ =3 γm − φ E (0) − ∂ ( Eµ ( E )) ∂E (¯ σ − σ ) µ ( E )(¯ σ − σ (0 , R )) = − ∂Q∂E | r =0 µ ( E )(¯ σ − σ (0)) > . It is a very heavy assumption to require ˜ θ < , and in some special cases of interest, this is notpossible. That is why we solve the ODE system starting from r = R. Therefore, we use a different approach. Instead of working from r = 0, we shall solve the ODEstarting from r = R. Lemma 4.8 (The bounds for E ) . There holds E ( R ) < E ( r ) < h ( r, R ) , r ∈ ( τ, R ) . Moreover, E ∈ C ( τ, R ] , E (cid:48) ( r ) < for r ∈ ( τ, R ) . Furthermore, if we define E ( τ ) = lim r → τ +0 E ( r ) ,then E ∈ C [ τ, R ] . Proof : Since the derivatives of E and h at r = R are given in (4.11) and (4.17), we deduce that d ( E − h ) dr (cid:12)(cid:12)(cid:12) r = R = (cid:20) Eµ ( E ) σ r Q E (1 , m, E ) − µ ( E )(1 − ¯ σ ) − Eµ ( E ) σ r Q E (1 , m, E ) (cid:21) (cid:12)(cid:12)(cid:12)(cid:12) r = R > . PARABOLIC-HYPERBOLIC SYSTEM MODELING 19
Since E ( R ) = h ( R ), we get E ( r ) < h ( r ) for 0 < R − r (cid:28) E ( r ) < h ( r ) for all r ∈ ( τ, R ). If there exists r ∈ ( τ, R ) such that E ( r ) = h ( r , R ) , E ( r ) < h ( r, R ) , r ∈ ( r , R ) , then E r ( r ) ≤ h r ( r ). However, h (cid:48) < r > E (cid:48) ( r ) = 0 by (4.16), this is a contradiction.Therefore, E ( r ) − h ( r ) < E (cid:48) ( r ) < r ∈ ( τ, R ), again by (4.16).From (4.4) or (4.12), we see that E is bounded and monotone decreasing in ( τ, R ]. Thus, the limitlim r → τ +0 E ( r ) exists and is positive, and is denoted by E ( τ ). Therefore, E ∈ C [ τ, R ] ∩ C ( τ, R ] . Now we set I ( r ) = r u ( r ) or I ( r ) (cid:44) (cid:90) rR µ ( E ( ρ ))( σ ( ρ ) − ¯ σ ) ρ dρ. As a consequence of Lemma 4.8, we find I ( r ) and its derivatives I (cid:48) ( r ) are continuous up to r = τ ,say, I ∈ C [ τ, R ]. Thus τ = τ ( R ) satisfies: τ = inf { τ (cid:48) ∈ (0 , R ) : I ( r ) < for r ∈ ( τ (cid:48) , R ) } . We have the following more precise results about values at r = τ . Lemma 4.9.
We have that(i) If τ ∈ (0 , R ) , then I ( τ ) = 0 , u ( τ ) = 0 . Moreover, E ( τ ) = h ( τ ) > , E ∈ C [ τ, R ] , E (cid:48) ( τ ) < , and E (cid:48) ( τ ) = E ( τ ) µ ( E ( τ )) σ r ( τ ) Q E ( σ ( τ ) , m, E ( τ )) − µ ( E ( τ ))( σ ( τ ) − ¯ σ ) . (4.19) (ii) If τ = 0 and I (0) = 0 , then u (0) = 0 , E (0) = h (0) , E ∈ C [0 , R ] and E (cid:48) (0) = 0 . Proof : The case τ > . Noting that ( τ, R ] is the maximal existence interval, we get u ( τ ) = 0 and I ( τ ) = 0. From I ( τ ) = I ( R ) = 0 and I (cid:48) ( r ) = µ ( E ( r ))( σ ( r ) − ¯ σ ) ρ , and σ r > , we see that σ ( r ) − ¯ σ = 0 has a root in [ τ, R ]. Moreover, the root is unique and in ( τ, R ), σ ( τ ) − ¯ σ < E ∈ C [ τ, R ] and u ( r ) = 1 r (cid:90) rR µ ( E )( σ ( ρ ) − ¯ σ ) ρ dρ = 1 r (cid:90) rτ µ ( E )( σ ( ρ ) − ¯ σ ) ρ dρ, we find u ∈ C [ τ, R ] and u ( τ ) = 0, u r ( τ ) = µ ( E ( τ ))( σ ( τ ) − ¯ σ ) < E ( τ ) (cid:54) = h ( τ ), then C = Q ( σ ( τ, R ) , m, E ( τ, R )) (cid:54) = 0. Setting C = u (cid:48) ( τ ) <
0, we deduce fromequation (4.16) that E (cid:48) = 1 r − τ C + o (1) C + o (1) near r = τ. Integrating the above equality, we see that E is unbounded near r = τ . This yields a contradiction.Hence, E ( τ ) = h ( τ, R ).From (4.4), we get θ = Q E ( σ ( τ ) , m, E ( τ )) u r ( τ ) = Q E ( σ ( τ ) , m, E ( τ )) µ ( E ( τ ))( σ ( τ ) − ¯ σ )= Q E ( σ ( τ ) , m, E ( τ )) − µ ( E ( τ ))( σ ( τ ) − ¯ σ ) µ ( E ( τ ))( σ ( τ ) − ¯ σ ) + 1 > . Now we apply Lemma 5.1 for equation (4.14a) to obtain that E ∈ C [ τ, R ] and E r ( τ ) = Q σ σ r u r − Q E | r = τ = E ( τ ) µ ( E ( τ )) σ r ( τ ) Q E ( σ ( τ ) , m, E ( τ )) − µ ( E )( σ ( τ ) − ¯ σ ) < . The case τ = 0 and I (0) = 0 . Noting that u ( r ) = 1 r (cid:90) rR µ ( E )( σ ( ρ ) − ¯ σ ) ρ dρ = 1 r (cid:90) r µ ( E )( σ ( ρ ) − ¯ σ ) ρ dρ, we get that u is differentiable at r = 0, and u (cid:48) (0) = µ ( E (0))( σ (0) − ¯ σ ) / < C (cid:54) = 0 where C = Q ( σ (0) , m, E (0)), then from equation (4.16), E (cid:48) = 1 r C + o (1) C + o (1) near r = 0 , where as before C = u (cid:48) (0) = µ ( E (0))( σ (0) − ¯ σ ) / <
0. This contradicts the boundedness of E andtherefore, C = Q ( σ (0) , m, E (0)) = 0, E (0) = h (0).From (4.4), we get θ = Q E (( σ (0) , m, E (0))) u r (0) = Q E (( σ (0) , m, E (0))) µ ( E (0))( σ (0) − ¯ σ ) /
3= 3 Q E (( σ (0) , m, E (0))) − µ ( E (0))( σ (0) − ¯ σ ) µ ( E (0))( σ (0) − ¯ σ ) + 1 > . Applying Lemma 5.1 for equation (4.14a), we obtain that E ∈ C [ τ, R ] and E (cid:48) (0) = 3 E (0) µ ( E (0)) σ r (0)3 Q E ( σ (0) , m ) − µ ( E (0))( σ (0) − ¯ σ ) = 0 . This completes the proof.
Remark 4.10.
The assumption (4.4) on Q , used in Lemma 4.9, only guarantee that the solution E obtained in Lemma 4.6 is differentiable at the minimal existence value r = τ . So the radiallysymmetric stationary solution obtained later in Theorem 4.17 satisfies E ∈ C [0 , . Continuous dependence with respect to the parameter R . In this subsection, we give some estimates of τ = τ ( R ) and the continuous dependence withrespect to R . Lemma 4.11. If < R (cid:28) , then τ ( R ) = 0 and I ( r, R ) < , ≤ r < R. Proof : Noting that 0 < ¯ σ <
1, andlim R → σ (0 , R ) − ¯ σ = lim R → R sinh( √ λR ) − ¯ σ = 1 − ¯ σ > , we see that for 0 < R (cid:28) σ ( r, R ) − ¯ σ ≥ σ (0 , R ) − ¯ σ > , r ∈ [0 , R ] . Observing that I ( R, R ) = 0 , I (cid:48) ( r, R ) = µ ( E )( σ ( r ) − ¯ σ ) r > , r ∈ [ τ ( R ) , R ) , we deduce that for 0 < R (cid:28) I ( r, R ) < , r ∈ [ τ ( R ) , R ) . The proof immediately follows from Lemma 4.9 (i).
Lemma 4.12. If R (cid:29) , then τ ( R ) > and I ( r, R ) < , τ ( R ) < r < R, I ( τ ( R ) , R ) = 0 . PARABOLIC-HYPERBOLIC SYSTEM MODELING 21
Proof : We take K such that e −√ λK = ¯ σ/ . Thenlim R →∞ σ ( R − K, R ) = lim R →∞ R sinh( √ λR ) sinh( √ λ ( R − K )) R − K = e −√ λK = 13 ¯ σ. Thus, for R (cid:29) , σ ( r, R ) − ¯ σ < σ ( R − K, R ) − ¯ σ <
12 ¯ σ, ≤ ρ ≤ R − K. From (4.12) and Lemma 4.8, we conclude E ( r, R ) is uniformly bounded, N < E ( r, R ) < N, r ∈ ( τ ( R ) , R ) , where N = E ( R, R ) = h ( R, R ) is independent of R . Thus, there exist two positive constants K < K (they are independent of R ) such that K < µ ( E ( r, R )) < K , r ∈ ( τ ( R ) , R ) . (4.20)Therefore, if τ ( R ) = 0 for some sufficiently large R , then0 ≤ (cid:90) R µ ( E )( σ − ¯ σ ) ρ dρ = (cid:90) R − K µ ( E )( σ − ¯ σ ) ρ dρ + (cid:90) RR − K µ ( E )( σ − ¯ σ ) ρ dρ< ( −
12 ¯ σ ) K ( R − K ) K ( R ) − ( R − K ) < , this is a contradiction. Thus, τ ( R ) > R . Lemma 4.13.
The function R ∈ (0 , ∞ ) (cid:55)→ τ ( R ) is upper semi-continuous, that is τ ( ¯ R ) ≥ lim sup R → ¯ R τ ( R ) for every ¯ R > . Proof : Let ¯
R > T ∈ ( τ ( ¯ R ) , ¯ R ) be any fixed constant. Then from Remark 4.7, wesee that the solution E ( r, R ) , u ( r, R ) of (4.14a), (4.14b) and (4.14c) exist at least for r ∈ [ T, R ) andfor R ∈ [ ¯ R − δ, ¯ R + δ ] for some δ >
0. Therefore, τ ( R ) satisfies τ ( R ) < T for every R ∈ [ ¯ R − δ, ¯ R + δ ] , and then lim sup R → ¯ R τ ( R ) ≤ T. Since T ∈ ( τ ( ¯ R ) , ¯ R ) is arbitrary, we finish the proof. Lemma 4.14.
The function τ ( R ) is lower semi-continuous, that is, for any ¯ R > , we have τ ( ¯ R ) ≤ lim inf R → ¯ R τ ( R ) . (4.21) Thus, the function R ∈ (0 , ∞ ) (cid:55)→ τ ( R ) is continuous. Proof : Noting that τ ( R ) ≥ R >
0, we know (4.21) holds when τ ( ¯ R ) = 0. Thus, we alwaysassume that ¯ τ = τ ( ¯ R ) >
0. Since I (¯ τ , ¯ R ) = 0 we find ψ (¯ τ , ¯ R ) < , where ψ ( r, R ) = σ ( r, R ) − ¯ σ . Bycontinuity there exists δ > ψ (¯ τ , ¯ R ) < ψ ( r, R ) < ψ (¯ τ , ¯ R ) < r ∈ [¯ τ − δ , ¯ τ + δ ] and R ∈ [ ¯ R − δ , ¯ R + δ ]. Suppose that (4.21) is false, then there exists δ > { R k } k ≥ with lim k →∞ R k =¯ R such that τ ( R k ) < ¯ τ − δ , k ≥ . (4.23)Without loss of generality, we may assume that 0 < δ < δ < ¯ τ and R k ∈ [ ¯ R − δ , ¯ R + δ ], k ≥
1. In order to yield a contradiction, we will show that I (¯ τ − δ , R k ) > k . By noting I (¯ τ , ¯ R ) = 0, we split I (¯ τ − δ , R k ) into three parts I (¯ τ − δ , R k ) =[ I (¯ τ − δ , R k ) − I (¯ τ + δ, R k )]+ [ I (¯ τ + δ, R k ) − I (¯ τ + δ, ¯ R )] + [ I (¯ τ + δ, ¯ R ) − I (¯ τ , ¯ R )] . From (4.3d), (4.20) and (4.22), we have I (¯ τ − δ , R k ) − I (¯ τ + δ, R k ) = − (cid:90) ¯ τ + δ ¯ τ − δ µ ( E ( r, R k )) ψ ( r, R k ) r dr ≥ − ψ (¯ τ , ¯ R )( δ + δ )(¯ τ − δ ) K > , (4.24) | I (¯ τ , ¯ R ) − I (¯ τ + δ, ¯ R ) | ≤ (cid:90) ¯ τ + δ ¯ τ µ ( E ( r, ¯ R )) | ψ ( r, ¯ R ) | r dr ≤ − ψ (¯ τ , ¯ R )(¯ τ + δ ) δK . (4.25)Note that E , u and I are continuous in ( r, R ) at ( r, R ) = (¯ τ + δ, ¯ R ), we derivelim k →∞ I (¯ τ + δ, R k ) − I (¯ τ + δ, ¯ R ) = 0 . (4.26)Combining (4.24), (4.25) and (4.26), we getlim inf k →∞ I (¯ τ − δ , R k ) ≥ − ψ (¯ τ , ¯ R )( δ + δ )(¯ τ − δ ) K + 2 ψ (¯ τ , ¯ R )(¯ τ + δ ) δK and by letting δ → k →∞ I (¯ τ − δ , R k ) ≥ − ψ (¯ τ , ¯ R ) δ (¯ τ − δ ) K > . This contradicts (4.23), By the definition of τ ( R k ) we conclude that I (¯ τ − δ , R k ) <
0. Thus (4.21)holds.Combining Lemma 4.13 and (4.21), we find that the function R ∈ (0 , ∞ ) (cid:55)→ τ ( R ) is continuous. Lemma 4.15.
The function I = I ( r, R ) is continuous in ( r, R ) for R > , τ ( R ) ≤ r ≤ R . Proof : Let ¯ R and ¯ r be fixed satisfying τ ( ¯ R ) ≤ ¯ r ≤ ¯ R . From (4.20), we obtain | I r ( r, R ) | ≤ K r , τ ( R ) ≤ r ≤ R, where the constant K is independent of r, R . For any (cid:15) >
0, we choose a constant ˆ δ such that16ˆ δK (¯ r + ˆ δ ) < (cid:15). Using Lemma 4.6 and Remark 4.7, we see I = I ( r, R ) is continuous at ( r, R ) = (¯ r + ˆ δ, ¯ R ) and thereis a constant δ ∈ (0 , ˆ δ ) such that for R ∈ [ ¯ R − δ, ¯ R + δ ],2 | I (¯ r + ˆ δ, R ) − I (¯ r + ˆ δ, ¯ R ) | < (cid:15). Therefore, for r ∈ [¯ r − δ, ¯ r + δ ] ∩ [ τ ( R ) , R ] and R ∈ [ ¯ R − δ, ¯ R + δ ], we have | I ( r, R ) − I (¯ r, ¯ R ) | ≤| I ( r, R ) − I (¯ r + ˆ δ, R ) | + | I (¯ r + ˆ δ, R ) − I (¯ r + ˆ δ, ¯ R ) | + | I (¯ r + ˆ δ, ¯ R ) − I (¯ r, ¯ R ) |≤ δK (¯ r + ˆ δ ) + | I (¯ r + ˆ δ, R ) − I (¯ r + ˆ δ, ¯ R ) | + 2ˆ δK (¯ r + ˆ δ ) ≤ δK (¯ r + ˆ δ ) + | I (¯ r + ˆ δ, R ) − I (¯ r + ˆ δ, ¯ R ) | < (cid:15) + 12 (cid:15) = (cid:15). PARABOLIC-HYPERBOLIC SYSTEM MODELING 23
It immediately follows that I is continuous at (¯ r, ¯ R ).4.4. Existence and Uniqueness for the stationary system.Lemma 4.16.
There exists a positive constant R ∗ such that τ ( R ∗ ) = 0 and I (0 , R ∗ ) = 0 . Proof : Let R ∗ be defined as R ∗ = inf { R > τ ( R ) > } . From Lemmas 4.11 and 4.12, we deduce that R ∗ is well-defined and is a positive finite number,and τ ( R ) = 0 for every 0 < R < R ∗ . By the continuity of R (cid:55)→ τ ( R ) (see Lemma 4.14), wederive τ ( R ∗ ) = 0. From the definition of R ∗ , there exists a sequence { R j } with R j > R ∗ andlim j →∞ R j = R ∗ such that τ ( R j ) > . From Lemma 4.9, E ( τ ( R j ) , R j ) = 0. It immediately followsfrom Lemmas 4.15 and 4.14 that E ( τ ( R ∗ ) , R ∗ ) = 0. This completes the proof.Combining Lemmas 4.16 and 4.9, we obtain the existence of radial stationary solution. Theorem 4.17.
The radially symmetric stationary problem (4.3) admits a solution.
The uniqueness.
In this subsection we establish the uniqueness of the radially symmetricstationary solution.We use a change of variables that transform the free boundary into a fixed boundary: s = rR , ˜ E ( s, R ) = E ( r, R ) , ˜ u ( s, R ) = r u ( r, R ) R . Then the initial problem (4.14a), (4.14b), (4.14c) of (
E, u ) is transformed to˜ u ˜ E s = s ( − γm ˜ E + φ ( ˜ E ) − ˜ Eµ ( ˜ E ) ˜ ψ ) , (4.27a)˜ u s = s µ ( ˜ E ) ˜ ψ, (4.27b)˜ u (1) = 0 , ˜ E (1) = ˜ h (1 , R ) , (4.27c)where ˜ h ( s ) = ˜ h ( s, R ) = h ( sR, R ), h ( r, R ) is solved from (4.10) and˜ ψ ( s, R ) = σ ( sR, R ) − ¯ σ = √ λR sinh( √ λR ) sinh( s √ λR ) s √ λR − ¯ σ. (4.28)Set ˜ τ ( R ) = τ ( R ) /R . Lemma 4.18.
For any R > R > u ( s, R ) > ˜ u ( s, R ) (4.29) for all > s ≥ max { ˜ τ ( R ) , ˜ τ ( R ) } . Moreover, ˜ τ ( R ) is an increasing function when ˜ τ ( R ) > . Proof : It is obvious that ˜ ψ (1 , R ), ˜ E (1 , R ), ˜ u s (1 , R ) are independent of R . One can easily computethe derivatives of ˜ ψ as follows ∂ R ˜ ψ ( s, R ) = √ λ cosh( √ λR ) cosh( s √ λR ) s [sinh( √ λR )] [ s tanh( √ λR ) − tanh( s √ λR )] < ,∂ s ˜ ψ (1 , R ) = √ λR coth( √ λR ) − > ,∂ sR ˜ ψ (1 , R ) = sinh(2 √ λR ) − √ λR ( √ λR ) √ λ > , (4.30) and calculate the derivatives of ˜ u and ˜ E as follows˜ u s (1 , R ) = µ ( ˜ E ) ˜ ψ | s =1 > , ˜ E s (1 , R ) = − ˜ Eµ ( ˜ E ) ˜ ψµ ( ˜ E ) ˜ ψ − Q E (cid:12)(cid:12)(cid:12)(cid:12) s =1 ∂ s ˜ ψ (1 , R ) < , ˜ E sR (1 , R ) = − ˜ Eµ ( ˜ E ) ˜ ψµ ( ˜ E ) ˜ ψ − Q E (cid:12)(cid:12)(cid:12)(cid:12) s =1 ∂ sR ˜ ψ (1 .R ) < . (4.31)Set ˜ E i ( s ) = ˜ E ( s, R i ), ˜ ψ i ( s ) = ˜ ψ ( s, R i ) , i = 1 ,
2. By (4.30) and (4.31), we get˜ E − ˜ E = 0 , [ ˜ E − ˜ E ] s < s = 1and hence [ ˜ E − ˜ E ]( s ) > < − s (cid:28)
1. Combining this with the fact that 0 < ˜ ψ < ˜ ψ ,0 < − s (cid:28) µ ( E ) is a decreasing positive function, we see that[˜ u − ˜ u ] s ( s ) = s [ µ ( ˜ E ( s )) ˜ ψ ( s ) − µ ( ˜ E ( s )) ˜ ψ ( s )] , = s [ µ ( ˜ E ( s )) − µ ( ˜ E ( s ))] ˜ ψ ( s ) + s µ ( ˜ E ( s ))[ ˜ ψ ( s ) − ˜ ψ ( s )] < u − ˜ u ]( s ) > < − s (cid:28)
1. Therefore,[ ˜ E − ˜ E ]( s ) > , [˜ u − ˜ u ]( s ) > < − s (cid:28) s, u − ˜ u ](¯ s ) >
0. Otherwise, ˜ u (¯ s ) = ˜ u (¯ s ) ≤
0. From (4.27), we obtain (˜ u ˜ E ) s =( − γm ˜ E + φ ( E )) s and ˜ u ( s ) ˜ E ( s ) = − (cid:90) s [ − γm ˜ E ( ξ ) + φ ( ˜ E ( ξ ))] ξ dξ. Thus, by using φ (cid:48) ≤
0, we derive0 ≥ ˜ u (¯ s )[ ˜ E − ˜ E ](¯ s )= − (cid:90) s { [ − γm ˜ E ( ξ ) + φ ( ˜ E ( ξ ))] − [ − γm ˜ E ( ξ ) + φ ( ˜ E ( ξ ))] } ξ dξ> . This yields a contradiction. In particular, ˜ u (¯ s ) < ˜ u (¯ s ) ≤ E − ˜ E ](¯ s ) >
0. Indeed, on the interval (¯ s,
1) we have ∂ s ˜ E ( s ) = Q ( ˜ ψ , m, ˜ E ) s ˜ u < Q ( ˜ ψ , m, ˜ E ) s ˜ u < Q ( ˜ ψ , m, ˜ E ) s ˜ u . It immediately follows from comparison principle of ordinary differential equation that ˜ E (¯ s ) > ˜ E (¯ s ).From the definition of ¯ s and the two assertions above, we conclude that ¯ s = max { ˜ τ ( R ) , ˜ τ ( R ) } and (4.29) hold. In particular, if ˜ τ ( R ) >
0, then by Lemma 4.9, ˜ u (˜ τ ( R ) , R ) = 0. Since ˜ u ( s, R ) < s ∈ (˜ τ ( R ) , R ), we get ˜ τ ( R ) > ˜ τ ( R ) > Theorem 4.19 (Uniqueness) . There exists exactly one
R ∈ (0 , ∞ ) such that τ ( R ) = 0 , I (0 , R ) = 0 ,i.e., (4.14) admit exactly one solution ( E , U , R ) . Proof : If there are two constants R > R > τ ( R i ) = 0 , I ( τ ( R i ) , R i ) = 0 , i = 1 , . PARABOLIC-HYPERBOLIC SYSTEM MODELING 25
Then ˜ τ ( R i ) = 0 and ˜ u (˜ τ ( R i ) , R i ) = 0. This contradicts Lemma 4.18 and hence the uniquenessfollows. 5. Numerical results and discussion
In this section, we investigate numerically the asymptotic stability of the system in the one-dimensional case. Note that when ECM is a constant, the stationary solution is stable for small µ and unstable for large µ . The uniqueness of our stationary solution indicates that our solutionshall be ”close” to the solution (corresponding to the constant ECM), when µ ( E ) is ”close” to aconstant. Hence the biological implication is that our stationary solution should be stable when µ ( E ) is small and unstable when µ ( E ) is large. However, it is a big challenge for us to confirm ourconjecture using mathematical analysis, since our system is very complex. Therefore, we performthe corresponding numerical simulations to confirm our expectation.Firstly, we solve the radially symmetric stationary equation (4.3), choose Q ( E, r ) = − γ αβ E + φ ( E ) − Eµ ( E )( √ λR sinh( √ λR ) ( sinh( √ λr ) √ λr − ˜ σ ) ,φ ( E ) = µ (1 − E ) , µ ( E ) = µ E , and R s = R to be determined. Then, we solve the system (2.9)and compare the long time behavior of E, σ, m with the steady state of
E, σ, m .To describe the long time behavior of our model, we simulate the time up to t = 200 , but formost of our simulations, the profile is already very close to steady state at t = 40, so we showmainly the dynamics for t ≤ µ ( E )is small, the stationary solution is stable, when µ ( E ) is big, the stationary solution is unstable.Therefore, we divided into three cases according to the the value of parameter µ . For each case, weinvestigate the dynamics of the density of ECM and concentration of the nutrient σ , but neglect theimpact of the density of MDE; as a matter of fact, we can prove that m ( x, t ) → αβ uniformaly as t →∞ .Note that m ( x, t ) satisfies ∂m ( x, t ) ∂t = D m (cid:52) m ( x, t ) + α − βm ( x, t ) , in Ω( t ) , (5.1a) ∂m ( x, t ) ∂n = 0 in ∂ Ω( t ) , (5.1b) m ( x,
0) = m ( x ) in Ω( t ) . (5.1c)We can compare m ( r, t ) with the solution ˆ m ( t ) of the ODE equation ∂ ˆ m ( t ) ∂t = α − β ˆ m ( t ) , ˆ m (0) = m , since ˆ m satisfies the same system (5.1) as m ( x, t ) and also satisfiesˆ m → αβ as t → ∞ , we deduce that m ( x, t ) → αβ uniformaly as t → ∞ , by a comparison theorem for parabolic equations [10]. Therefore, it sufficed to consider the specialcase m = αβ from now on. Table 1. parameter ranges in the modelParameters Description Range Reference c Nutrient diffusion coefficient 10 − − − [4, 24] λ proliferation of nutrient 0.05-2 [4] D m MDE diffusion coefficient 10 − −
10 [4, 24] µ mobility coefficient 0.9-1.45 [4, 24] µ proliferation of ECM coefficient 0.15-2.5 [4] γ Rate of degradation of ECM 1-20 [4, 24] α Production of MDEs 0.01-5 [24] β Decay of MDE 10 − −
10 [24]5.1.
Case I: µ = 0 . . Figure 1.
Figure 1: One dimensional numerical results for case 1 when µ = 0 .
5. Results are snapshotsof the stationary system for
E, u and σ at γ = 10; α = 0 . β = 1; λ = 2; ˜ σ = 0 . µ = 0 . R =1 . m = αβ . The horizontal axis r indicates the spatial position, and the vertical axis indicatesthe density of ECM, the velocity or the concentration of nutrient listed in the legend.In order to show the profile of the evolution of the concentration of the nutrient σ , the densityof ECM and MDE, firstly, we chose the steady state to be the initial conditions for all the system,which is a perfect case. PARABOLIC-HYPERBOLIC SYSTEM MODELING 27
Figure 2.
Next we consider the initial condition to be a small perturbation of the steady state.
Figure 3.
Figure 3: In order to verify the stability of the steady state solution, we simulate a large numberof initial conditions and observe the long time behavior. Under these conditions, our simulation shows, when t −→ ∞ σ ( r, t ) −→ the steady state: sinh √ λrr R sinh √ λR ,E ( r, t ) −→ the steady state: E s ,m ( r, t ) −→ the steady state: αβ . These plots are for the case µ = 0 . Case II: µ = 3 . . Figure 4.
Figure 4: One dimensional numerical results for case 2 when µ = 3 .
0. Results are snapshotsof the stationary system for
E, u and σ at γ = 10; α = 0 . β = 1; λ = 2; ˜ σ = 0 . µ = 0 . R =1 . m = αβ = 0 . . The horizontal axis r indicates the spatial position, and the vertical axisindicates the density of ECM, the velocity or the concentration of nutrient listed in the legend.
PARABOLIC-HYPERBOLIC SYSTEM MODELING 29
Figure 5.
Figure 5: In order to show the profile of the evolution of the concentration of the nutrient σ , thedensity of ECM and MDE, firstly, we chose the steady state to be the initial conditions for all thesystem, which is a perfect case.Next we consider the initial condition to be a small perturbation of the steady state. Figure 6.
Figure 6: In order to verify the stability of the steady state, we simulate a large number of initialconditions and observe the long time behavior. Under these conditions, our simulation shows, when t −→ ∞ σ ( r, t ) −→ the steady state: sinh √ λrr R sinh √ λR ,E ( r, t ) −→ the steady state: E s ,m ( r, t ) −→ the steady state: αβ . These plots are for the case µ = 3 . Case III: µ = 10 . We now consider the parameter µ = 10, and obtain the critical R ≈ . σ , thedensity of ECM and MDE, we still consider the initial condition to be a small perturbation of thesteady state. We found the steady state of E is about to stay at 0 . − . µ = 10, however,the E ( r, t ) goes to a value range from 0 . − . Figure 7.
Figure 7: One dimensional numerical results for case 3 when µ = 10 .
0. Results are snapshotsof the stationary system for
E, u and σ at γ = 10; α = 0 . β = 1; λ = 2; ˜ σ = 0 . µ = 0 . R =1 . m = αβ = 0 . . The horizontal axis r indicates the spatial position, and the vertical axisindicates the density of ECM, the velocity or the concentration of nutrient listed in the legend.
Appendix 1: Proof of Theorem 4.1
For simplicity, we may assume that x = 0 and f (0 ,
0) = 0 , g (0 ,
0) = 1 , θ = f x (0 , < . (5.2)To prove Theorem 4.1, we shall work on the solution for positive time t >
0. The solution fornegative time t <
Step 1.
Existence.In order to show the existence of solution of (4.1), we shall approximate the singular equationwith non-singular ones. The approximated solution x (cid:15) , (cid:15) > PARABOLIC-HYPERBOLIC SYSTEM MODELING 31 dx (cid:15) dt = f ( x (cid:15) , t ) (cid:90) (cid:15) g (0 , ρ ) dρ + (cid:90) t(cid:15) g ( x (cid:15) , s ) ds , r > (cid:15),x (cid:15) ( t ) = 0 , ≤ t ≤ (cid:15). (5.3)For (cid:15) sufficiently small, (cid:82) (cid:15) g (0 , ρ ) dρ > dx (cid:15) /dr doesnot vanish at r = (cid:15) . For such an (cid:15) , the equation (5.3) does not have singularity. If one denotesthat denominator by y = y ( t ), then x = x ( t ) , y = y ( t ) satisfies ordinary differential equations. Bythe standard ODE theory, the solution x (cid:15) of (5.3) exists and is unique, and the maximal existenceinterval is denoted by [0 , T (cid:15) ) for 0 < (cid:15) (cid:28)
1. Moreover, x (cid:15) satisfies x (cid:15) ( t ) = (cid:90) t f ( x (cid:15) ( τ )) , τ ) (cid:82) τ g ( x (cid:15) ( s ) , s ) ds dτ, t > (cid:15). (5.4)We first estimate the lower bound and the upper bound of x (cid:15) .Let η > η < − θ and η <
1. By the continuity, there exists a δ = δ ( η ) > | g ( x, t ) − | < η, | f x ( x, t ) − θ | < η, | f ( x, t ) | < η (5.5)for all ( x, t ) ∈ ¯ D δ , where ¯ D δ = { ( x, t ) : | x | ≤ δ, | t | ≤ δ } . Let M > | f x ( x, t ) | < M, | f t ( x, t ) | < M for ( x, t ) ∈ ¯ D δ . Now we take two positive constants κ and T such that κ > M − θ − η , T = min { δ, δκ } . (5.6)We claim that x (cid:15) exists for t ∈ [0 , T ] for every (cid:15) ∈ (0 , T ), − κt < x (cid:15) ( t ) < κt, ≤ t ≤ T , (5.7) | x (cid:48) (cid:15) ( t ) | ≤ M κ + M − η , (cid:15) < t ≤ T . (5.8)Indeed, we have for 0 < t ≤ T , κ (1 − η ) t − f ( κt, t ) > κ (1 − η ) t − ( θ + η ) κt − M t = [(1 − θ − η ) κ − M ] t > , − κ (1 − η ) t − f ( − κt, t ) < − κ (1 − η ) t + ( θ + η ) κt + M t = − [(1 − θ − η ) κ − M ] t < . Thus, x ( t ) = κt and x ( t ) = − κt satisfy x (cid:48) ( t ) > f ( x ( t ) , t )(1 − η ) t , x (cid:48) ( t ) < f ( x ( t ) , t )(1 − η ) t , < t ≤ T . Now we proceed to prove (5.7). If there is a first ¯ t ∈ ( (cid:15), T ] ∩ ( (cid:15), T (cid:15) ) such that x (cid:15) ( t ) < x ( t ) , t ∈ ( (cid:15), ¯ t ) and x (cid:15) (¯ t ) = x (¯ t ) , then, by the definition of derivatives, x (cid:48) (cid:15) (¯ t ) ≥ x (cid:48) (¯ t ) = κ >
0, and f ( x (cid:15) (¯ t ) , ¯ t ) > x (cid:15) . However, at t = ¯ t , x (cid:48) (cid:15) (¯ t ) = f ( x (cid:15) (¯ t ) , ¯ t ) (cid:82) ¯ t g ( x ( s ) , s ) ds < f ( x (¯ t ) , ¯ t )(1 − η )¯ t < x (cid:48) (¯ t ) . This yields a contradiction. Thus, x (cid:15) ( t ) < x ( t ), t ∈ ( (cid:15), T ] ∩ ( (cid:15), T (cid:15) ). Similarly, we can prove x (cid:15) > x ( t ), t ∈ ( (cid:15), T ] ∩ ( (cid:15), T (cid:15) ). Since T (cid:15) is the maximal existence time of x (cid:15) , we see T (cid:15) > T . Therefore, (5.7)and then (5.8) are established. Now we show the existence of solution of (4.1). Indeed, from (5.7) and (5.8), we see that x (cid:15) is uniformly bounded in the Lipschitz function space Lip[0 , T ]. By passing to a subsequence ifnecessary, we may assume that x (cid:15) → x in C [0 , T ]as (cid:15) ↓ x ∈ Lip[0 , T ]. From (5.4), we obtain x ∈ Lip[0 , T ] satisfies x ( t ) = (cid:90) t f ( x ( τ ) , τ ) (cid:82) τ g ( x ( s ) , s ) ds dτ. Thus, x ∈ Lip[0 , T ] ∩ C (0 , T ] is a solution to (4.1). Moreover, x is of class C [0 , T ] and x (cid:48) (0) = k = γ/ (1 − θ ), γ = f t (0 , Step 2.
Uniqueness when θ ∈ ( − , x ( t ) , y ( t ) are two C solutions of (4.1) and they are well-defined at a commoninterval [0 , T ]. Set L = max { max t ∈ [0 ,T ] | x (cid:48) ( t ) | , max t ∈ [0 ,T ] | y (cid:48) ( t ) |} . (5.9)Let us fix η ∈ (0 ,
1) satisfying ˜ θ = | θ | + η − η ∈ (0 , . (5.10)By the continuity, there exists a δ = δ ( η ) > M = M δ > | f x ( x, t ) | < M, | f t ( x, t ) | < M, | g ( x, t ) − g ( y, t ) | ≤ M | x − y | (5.11)for all ( x, t ) and ( y, t ) in ¯ D δ . Thus, z ( t ) = x ( t ) − y ( t ) satisfies dz ( t ) dt = f ( x ( t ) , t ) − f ( y ( t ) , t ) (cid:82) t g ( x ( s ) , s ) ds − f ( y ( t ) , t ) (cid:82) t [ g ( x ( s ) , s ) − g ( y ( s ) , s )] ds (cid:82) t g ( x ( s ) , s ) ds (cid:82) t g ( y ( s ) , s ) ds , it follows from (5.5) and (5.11) that | z (cid:48) ( t ) | ≤ ( | θ | + η ) | z ( t ) | (1 − η ) t + M ( L + 1) t · tM (1 − η ) t max s ∈ [0 ,t ] | z ( s ) | , for 0 < t ≤ T = min { δ, δL } . If we set ˜ M = M ( L + 1)(1 − η ) − and Z ( t ) = max s ∈ [0 ,t ] | z ( s ) | , then | z (cid:48) ( s ) | ≤ ˜ θ | z ( s ) | s + ˜ M Z ( t ) , (5.12)for all s ∈ (0 , t ] and all t ∈ (0 , T ].We claim that 0 ≤ | z ( s ) | ≤ ˜ M s − ˜ θ Z ( t ) , s ∈ (0 , t ] , t ∈ (0 , T ] . (5.13)Indeed, from (5.9) and (5.12), one easily obtain by induction that | z (cid:48) ( s ) | ≤ a k , | z ( s ) | ≤ a k t, s ∈ [0 , t ] . Here a = 2 L and a k +1 = ˜ θa k + ˜ M Z ( t ). (5.13) follows by the following factlim k →∞ a k = ˜ M − ˜ θ Z ( t ) . Take t = t ∗ ∈ (0 , T ] such that ˜ M t ∗ < − ˜ θ , then0 ≤ Z ( t ) ≤ ˜ M t ∗ − ˜ θ Z ( t ) , t ∈ [0 , t ∗ ] . PARABOLIC-HYPERBOLIC SYSTEM MODELING 33
It immediately follows that Z ( t ∗ ) = 0 and hence z ( t ) = 0, 0 ≤ t ≤ t ∗ . This shows the uniqueness inthe case of θ ∈ ( − , Step 3.
Uniqueness when θ ∈ ( −∞ , x ( t ) , y ( t ) are two solutions of (4.1) and they are well-defined at a common interval[0 , T ]. Then z ( t ) = x ( t ) − y ( t ) satisfies( z ( t ) t − θ ) (cid:48) = t − θ [ z (cid:48) ( t ) − θz ( t ) t ] =: t − θ A ( t ) , where A ( t ) = (cid:32) f ( x ( t ) , t ) (cid:82) t g ( x ( s ) , s ) ds − θx ( t ) t (cid:33) − (cid:32) f ( y ( t ) , t ) (cid:82) t g ( y ( s ) , s ) ds − θy ( t ) t (cid:33) = [ f ( x ( t ) , t ) − θx ( t )] − [ f ( y ( t ) , t ) − θy ( t )] (cid:82) t g ( x ( s ) , s ) ds + θ [ x ( t ) − y ( t )] (cid:82) t [1 − g ( x ( s ) , s )] dst (cid:82) t g ( x ( s ) , s ) ds − f ( y ( t ) , t ) (cid:82) t [ g ( x ( s ) , s ) − g ( y ( s ) , s )] ds (cid:82) t g ( x ( s ) , s ) ds (cid:82) t g ( y ( s ) , s ) ds . Let δ ∈ (0 , T ) and M > | g ( x, t ) − g ( y, t ) | ≤ M | x − y | (5.14)for all ( x, t ) and ( y, t ) in ¯ D δ . Let us fix η > δ ∈ (0 , δ ) > α ∈ (0 ,
1) with α = 1 | θ | (cid:18) η − η + | θ | η (1 − η ) + M η (1 − η ) (cid:19) (5.15)and that (5.5) holds. Now let t ∗ ∈ (0 , δ ) such that | x ( t ) | < δ, | y ( t ) | < δ, t ∈ [0 , t ∗ ] . Using (5.5), (5.14) and (5.15), we see that for t ∈ (0 , t ∗ ], | A ( t ) | ≤ η | z ( t ) | (1 − η ) t + | θ || z ( t ) | · tηt (1 − η ) + η · M Z ( t ) t (1 − η ) t ≤ α | θ | Z ( t ) t , where Z ( t ) = max s ∈ [0 ,t ] | z ( s ) | is a nondecreasing, nonnegative function. Thus, | ( z ( t ) t − θ ) (cid:48) | ≤ − θαt − θ − Z ( t ) , t ∈ (0 , t ∗ ] . Integrating the inequality above, we get | z ( t ) t − θ | ≤ − θαZ ( t ) (cid:90) t s − θ − ds = αZ ( t ) t − θ . It immediately follows that 0 ≤ Z ( t ) ≤ αZ ( t ), t ∈ [0 , t ∗ ]. Hence Z ( t ) = 0, x ( t ) − y ( t ) = z ( t ) = 0, t ∈ [0 , t ∗ ]. This shows the uniqueness in the case of θ ∈ ( −∞ , Appendix 2: Proof of Theorem 4.4:
Without loss of generality, we may assume that g ( ϕ ( µ ) , , µ ) ≡ f x ( ϕ ( µ ) , , µ ) = θ ( µ ) . Step 1.
We prove that there exists two constants T > (cid:15) > µ ∈ B µ ,(cid:15) , the solution x = x ( t, µ ) of (4.2) exists in [0 , T ];(ii) The function ( t, µ ) (cid:55)→ x ( t, µ ) is continuous in [0 , T ] × B µ ,(cid:15) . Since θ = θ (0) < θ ( µ ) is continuous for µ in a neighborhood of µ , there exists a δ such that θ ( µ ) < | µ − µ | ≤ δ . Let us fix η ∈ (0 ,
1) satisfyinginf { − θ ( µ ) − η : µ ∈ B µ ,δ } > . For such an η , there exists a δ ∈ (0 , δ ) and M > | f x ( x, t, µ ) − θ ( µ ) | < η, | g ( x, t, µ ) − | < η, | f x ( x, t, µ ) | < M, | f x ( x, t, µ ) | < M for all ( x, t, µ ) satisfying | x − ϕ ( µ ) | < δ , | t | < δ , µ ∈ B µ ,δ . Now we fix κ and T > κ > max { M − θ ( µ ) − η : µ ∈ B µ ,δ } and T = min { δ , δ κ } . From Theorem 4.1, the solution x ( t, µ ) of (4.2) exists and is unique, the (right) maximal existenceinterval is denoted by [0 , T µ ). From the proof of the existence results (see step 1 in the proof ofTheorem 4.1), we see that T µ > T , | x ( t, µ ) − ϕ ( µ ) | < κt, (cid:12)(cid:12)(cid:12)(cid:12) ∂x ( t, µ ) ∂t (cid:12)(cid:12)(cid:12)(cid:12) ≤ M κ + M − η for all t ∈ [0 , T ] and µ ∈ B µ ,δ . Thus { x ( · , µ ) : µ ∈ B µ ,δ } is a bounded in C [0 , T ]. By the Ascoli-Arzel´a theorem, { x ( · , µ ) : µ ∈ B µ ,δ } is a pre-compact subset in C [0 , T ]. Then for any ¯ µ ∈ B µ ,δ and convergent sequence { µ k } with limit ¯ µ , there exists a further subsequence { µ k l } such that thecorresponding subsequence of function { x ( · , µ k l ) } is a convergent subsequence in C [0 , T ], the limitis denoted ¯ x , which belongs to Lip[0 , T ]. Note that the solution x ( t, µ ) satisfies x ( t, µ ) = ϕ ( µ ) + (cid:90) t f ( x ( t, µ )) , t, µ ) (cid:82) t g ( x ( s, µ ) , s, µ ) ds , t > . (5.16)We know the limit function ¯ x ∈ Lip[0 , T ] also satisfies (5.16) with µ = ¯ µ . Hence ¯ x is the uniquesolution of (4.2) with µ = ¯ µ , the limit function ¯ x = x ( · , ¯ µ ) depends only on ¯ µ , which is independentof the choice of subsequence of { µ k } . Thus, the full sequence { x ( · , µ k ) } converges to x ( · , ¯ µ ) in C [0 , T ]. Therefore, ( t, µ ) (cid:55)→ x ( t, µ ) is continuous in [0 , T ] × B µ ,δ . Step 2.
We complete the proof. We rewrite our equation in the form dxdt = f ( x,t,µ ) y , dydt = g ( x, t, µ ) ,x (0) = 0 , y (0) = 0 . (5.17)Theorem 4.1 and step 1 suggest that the C solution ( x ( t, µ ) , y ( t, µ )) of (5.17) exists and is unique ina time interval [0 , T ] for µ ∈ B µ ,δ , and the solution is continuous in ( t, µ ) for t ∈ [0 , T ], µ ∈ B µ ,δ , y ( t, µ ) > (1 − η ) t , ( t, µ ) ∈ (0 , T ] × B µ ,δ . Moreover, when µ = µ , the solution ( x ( t, µ ) , y ( t, µ ))exists in a bounded closed interval [0 , b ].In order to show the continuous dependence in µ beyond time τ , we consider the following initialproblem d ˆ xdt = f (ˆ x,t,µ ) y , d ˆ ydt = g (ˆ x, t, µ ) , ˆ x ( T ) = ξ, ˆ y ( T ) = ζ. (5.18)Set ξ = x ( T , µ ), ζ = y ( T , µ ). The solution of (5.18) exists and is denoted by ˆ x = ˆ x ( t, ξ, ζ, µ ) , ˆ y =ˆ y ( t, ξ, ζ, µ ) for at least ( ξ, ζ, µ ) in a neighborhood of ( ξ , ζ , µ ). By the classical ODE theory,solutions of (5.18) and (5.17) are the same when ξ = x ( T , µ ), ξ = y ( T , µ ). In particular, when( ξ, ζ, µ ) = ( ξ , ζ , µ ), the solution (ˆ x, ˆ y ) exists in a bounded closed interval [ T , b ]. From the classicalODE theory, there exists a δ ∈ (0 , δ ) such that (i) the solution ˆ x = ˆ x ( t, ξ, ζ, µ ) , ˆ y = ˆ y ( t, ξ, ζ, µ ) of(5.18) exist in [ T , b ] for every | ξ − ξ | < δ , | ζ − ζ | < δ and µ ∈ B µ ,δ ; (ii) ˆ x = ˆ x ( t, ξ, ζ, µ ) , ˆ y = PARABOLIC-HYPERBOLIC SYSTEM MODELING 35 ˆ y ( t, ξ, ζ, µ ) are continuous functions in ( t, ξ, ζ, µ ) ∈ [ T , b ] × [ ξ − δ , ξ + δ ] × [ ζ − δ , ζ + δ ] × B µ ,δ .By step 1, we know there exists a δ ∈ (0 , δ ) such that | x ( T , µ ) − x ( T , µ ) | < δ , | y ( T , µ ) − y ( T , µ ) | < δ , µ ∈ B µ ,δ . Thus, x ( t, µ ) = ˆ x ( t, x ( T , µ ) , y ( T , µ ) , µ ) and y ( t, µ ) = ˆ y ( t, x ( T , µ ) , y ( T , µ ) , µ ) exist in the timeinterval [ T , b ] for every µ ∈ B µ ,δ , and these functions are continuous in ( t, µ ) for t ∈ [ T , b ] and µ ∈ B µ ,δ . Therefore, the solution x ( t, µ ) of (4.2) exists in [0 , b ] for µ ∈ B µ ,δ , and x ( t, µ ) iscontinuous in ( t, µ ) ∈ [0 , b ] × B µ ,δ .Similarly, the solution x ( t, µ ) of (4.2) exists in [ a,
0] for µ ∈ B µ ,δ , and x ( t, µ ) is continuous in( t, µ ) ∈ [ a, × B µ ,δ for some δ . Theorem 4.4 is proved by taking δ = min { δ , δ } . Appendix 3: Fixing an error in [12]There is an important reference [12], which is also a singular equation but without the nonlocalintegral term. However, we found an error in the proof of Lemma 9.3 in [12] and correct it asfollows:
Lemma 5.1 ([12]) . Let φ ∈ C [0 , a ] ∩ C (0 , a ] be a solution of v ( r ) dφdr = F ( φ, r ) , where v ∈ C [0 , a ] and F ∈ C in an open set containing the image { ( φ ( r ) , r ) : 0 ≤ r ≤ a } andsatisfies v (0) = 0 , F ( φ (0) ,
0) = 0 and F φ ( φ (0) , (cid:54) = v (cid:48) (0) , v (cid:48) (0) (cid:54) = 0 . (5.19) Then we have that φ ∈ C [0 , a ] and lim r ↓ dφdr ( r ) = dφdr (0) = F r ( φ (0) , v (cid:48) (0) − F φ ( φ (0) , . (5.20)This is Lemma 9.3 in [12].Consider the following equation (cid:40) dφdr = F ( φ,r ) r , r > ,φ (0) = 0 , (5.21)where F is a smooth function with F (0 ,
0) = 0.
Example 1. F ( φ, r ) = αφ with α ∈ (0 , φ ( r ) = cr α , r ≥ . Example 2. F ( φ, r ) = φ . Then (5.21) has exactly two kind of continuous solution: φ ( r ) ≡ φ ( r ) = 1ln( c/r ) , < r < c These two examples satisfy the condition (5.19), while the continuous solutions of (5.21) is notunique and is not differentiable at r = 0 except the trial solution φ ( r ) ≡ F φ ( φ (0) , /v (cid:48) (0) ∈ [0 , v (cid:48) (0) (cid:54) = 0. Proposition 5.2.
Let the assumptions in Lemma 5.1 hold. Set γ = F φ ( φ (0) , /v (cid:48) (0) . Then(i) If γ (cid:54)∈ [0 , , then φ ∈ C [0 , a ] and (5.20) holds.(ii) If γ ∈ [0 , and φ ∈ Lip [0 , a ] , then φ ∈ C [0 , a ] and (5.20) holds. Proof : As in the appendix in [12], we rewrite our equation in the linear form as follows dφdr = α ( r )( φ − φ (0)) + α ( r ) , where α ( r ) = F ( φ ( r ) , r ) − F ( φ (0) , r )( φ ( r ) − φ (0)) v ( r ) = F φ ( φ (0) , v (cid:48) (0) r (1 + o (1)) ,α ( r ) = F ( φ (0) , r ) − F ( φ (0) , v ( r ) = F r ( φ (0) , v (cid:48) (0) (1 + o (1)) , we then have ddr { ( φ ( r ) − φ (0)) G ( r ) } = α ( r ) G ( r ) (5.22)with G ( r ) = exp (cid:18)(cid:90) ar α ( s ) ds (cid:19) , G (cid:48) ( r ) = − α ( r ) G ( r ) . Case 1. F φ ( φ (0) , /v (cid:48) (0) > , v (cid:48) (0) (cid:54) = 0. This case is proved in the appendix of [12]. Case 2. F φ ( φ (0) , /v (cid:48) (0) < , v (cid:48) (0) (cid:54) = 0. In this case, there exist (cid:15) > r (cid:15) ∈ (0 , a ) such that rα ( r ) < − (cid:15) for 0 < r < r (cid:15) . Then 0 < G ( r ) < C (cid:15) r (cid:15) , r ∈ (0 , a ]for some C (cid:15) >
0. We integrate (5.22) over [0 , r ] to obtain( φ ( r ) − φ (0)) G ( r ) = (cid:90) r α ( η ) G ( η ) dη. (5.23)By L’Hˆospital’s rule, lim r ↓ φ ( r ) − φ (0) r = lim r ↓ (cid:82) r α ( η ) G ( η ) dηrG ( r )= lim r ↓ α ( r ) G ( r )[1 − rα ( r )] G ( r )= lim r ↓ α ( r )1 − rα ( r ) = F r ( φ (0) , v r (0) − F φ ( φ (0) ,
0) (5.24)and lim r ↓ dφ ( r ) dr = lim r ↓ (cid:20) rα ( r ) φ ( r ) − φ (0) r + α ( r ) (cid:21) = F φ ( φ (0) , v (cid:48) (0) F r ( φ (0) , v r (0) − F φ ( φ (0) ,
0) + F r ( φ (0) , v (cid:48) (0)= F r ( φ (0) , v r (0) − F φ ( φ (0) , . (5.25)Therefore, φ ∈ C [0 , a ] and (5.20) holds. Case 3. γ ∈ ( −∞ ,
1) and φ ∈ Lip [0 , a ]. In this case, rα ( r ) < − (cid:15) for 0 < r < r (cid:15) for some (cid:15) > r (cid:15) ∈ (0 , a ). Then 0 < G ( r ) < C (cid:15) r (cid:15) − , r ∈ (0 , a ]for some C (cid:15) >
0. Thus the right hand side of (5.22) is integrable near r = 0. Since φ ∈ Lip [0 , a ],we know | φ ( r ) − φ (0) | ≤ Cr , r ∈ [0 ,
1] for some C . We knowlim r → ( φ ( r ) − φ (0)) G ( r ) = 0 . We now integrate (5.22) over [0 , r ] to obtain (5.23). Thus, (5.24) and (5.25) is also obtained.Therefore, φ ∈ C [0 , a ] and (5.20) holds. Remark 5.3. If γ ∈ (0 , , then φ ∈ C γ − (cid:15) [0 , a ] for any small (cid:15) > ; If F φ ( φ (0) , (cid:54) = 0 and v ( r ) = cr β (1 + o (1)) for some c (cid:54) = 0 , β > , then φ has a derivative at r = 0 . PARABOLIC-HYPERBOLIC SYSTEM MODELING 37
Acknowledgments:
The authors are very grateful to the professor Yuan Lou for his helpfulcomments and Wenrui Hao for his useful suggestions. Rui Li is sponsored by the China ScholarshipCouncil. Rui Li also would like to thank the Department of Applied Computational Mathematicsand Statistics of the University of Notre Dame for its hospitality when she was a visiting student.
References [1] A.R.A. Anderson, M.A.J. Chaplain, E.L. Newman, R.J.C. Steele, A.M. Thompson, Mathematical modelling oftumour invasion and metastasis, J. Theor. Med. 2(2000), 129-154.[2] V. Andasari, A. Gerisch, G. Lolas, A. South and M.A.J. Chaplain, Mathematical modeling of cancer cell invasionof tissue: Biological insight from mathematical analysis and computational simulation, J. Math. Biol. 63(2011),141-172.[3] Frantz, C., Stewart, K.M., Weaver,V.M., The extracellular matrix at a glance. J. Cell Sci. 123(2010), 4195-4200.[4] M.A.Chaplain, G.Lolas, Modelling cell movement in anisotropic and heterogeneous tissue: dynamic heterogene-ity, Networks and Media 1(2016), 399-439.[5] S. Cui, A. Friedman, A hyperbolic free boundary problem modeling tumor growth, Interfaces Free Bound.5(2003), 159-181.[6] E. DiBenedetto, C.M. Elliott, A. Friedman, The free boundary of a flow in a porous body heated from itsboundary. Nonlinear Anal. 10(1986), 879-900.[7] P. Domschke, D. Trucu, A. Gerisch, M.A.J. Chaplain, Mathematical modelling of cancer invasion: Implicationsof cell adhesion variability for tumour infiltrative growth patterns, J. Theoret. Biol. 361(2014), 41-60.[8] H. Enderling, A.R.A. Anderson, M.A.J. Chaplain, G.W.A. Rowe, Visualisation of the numerical solution ofpartial differential equation systems in three space dimension and its importance for mathematical models inbiology, Math. Biosci. Eng. 3(2006), 571-582.[9] H. Enderling, A.R.A. Anderson, Mathematical Modeling of Tumor Growth and Treatment. Curr. Pharm. Des.20(2014), 1-7.[10] A. Friedman, Partial differential equations of parabolic type, Prentice-Hall, Inc., Englewood Cliffs, N.,J., 1964.[11] A. Friedman, Free boundary problems arising in tumor models. Atti Accad. Naz. Lincei Cl. Sci. Fis. Mat. Natur.Rend. Lincei, Mat. Appl. 15(2004), 161-168.[12] A. Friedman, B. Hu, The role of oxygen in tissue maintenance: mathematical modeling and qualitative analysis,Math. Models Methods Appl. Sci. 18(2008), 1409-1441.[13] A. Friedman, B. Hu, A Stefan problem for a protocell model. SIAM J. Math. Anal. 30 (1999), 912-926.[14] A. Friedman, K.-Y. Lam, Analysis of a free-boundary tumor model with angiougenesis. J. Differential Equations259(2015), 7636-7661.[15] A.Fredman, F. Reitich, Analysis of a mathematical model for the growth of tumors, J. Math.Biol. 38(1999),262-284.[16] A. Gerisch, M.A.J. Chaplain, Mathematical modelling of cancer cell invasion of tissue: Local and non-localmodels and the effect of adhesion, J. Theoret. Biol. 250(2008), 684-704.[17] H.P. Greenspan, Models for the growth of a solid tumor by diffusion, Stud. Appl. Math. 52(1972), 317-340.[18] G.M. Lieberman, Second Order Parabolic Differential Equations, World Scientific Publishing Co., Inc., RiverEdge, NJ, 1996.[19] J.S. Lowengrub, H.B. Frieboes, F. Jin, Y.L. Chuang, X. Li, P. Macklin, S.M. Wise, V. Cristini, Nonlinearmodelling of cancer: bridging the gap between cells and tumours, Nonlinearity 23(2010), R1-R9.[20] O.A. Ladyˇzenskaja, V.A. Solonnikov, N.N. Ura´lceva, Linear and Quasi-linear Equations of Parabolic Type,American Mathematical Society, Providence, 1968.[21] L. Peng, D. Trucu, A. Thompson, P. Lin, M.A.J. Chaplain, A multiscale mathematical model of tumour invasivegrowth, Bull. Math. Biol. 79(2017), 389-429.[22] Thomas R. Cox, Janine T. Erler, Remodeling and homeostasis of the extracellular matrix: implications forfibrotic diseases and cancer, Dis. Model Mech. 4(2011), 165-178.[23] Stetler Stevenson, W. G., Aznavoorian, S. Liotta L. A., Tumor cell interactions with the extracellular matrixduring invasion and metastasis, Annu. Rev. Cell Biol. 9(1993), 541-573.[24] Z. Szymanska, J. Urbanski, A. Marciniak-Czochra, Mathematical modelling of the influence of heat shock proteinson cancer invasion of tissue, J. Math. Biol. 58(2009), 819-844.[25] D. Trucu, P. Lin, M.A.J. Chaplain, Y. Wang, A multiscale moving boundary model arising in cancer invasion,Multiscale Model. Simul. 11(2013), 309-335.[26] Y. Tao, A free boundary problem modeling the cell cycle and cell movement in multicellular tumor spheroids,J. Differential Equations 247(2009), 49-68. [27] Yen T. Nguyen Edalgo, Ashlee N. Ford Versypt, Mathematical modeling of metastatic cancer migration througha remodeling extracellular matrix, Processes 6(2018),1-16.[28] J. Wu, F. Zhou, Bifurcation analysis of a free boundary problem modelling tumour growth under the action ofinhibitors, Nonlinearity 25(2012), 2971-2991.[29] L.B. Weiswald, D.Bellet, V. Dangles-Marie, Spherical cancer models in tumor biology, Neoplasia 17(2015), 1-15.[30] M. Wu, H.B. Frieboes, M.A.J. Chaplain, S. McDougall, V. Cristini, J.S. Lowengrub, The effect of interstitialpressure on therapeutic agent transport: coupling with the tumor blood and lymphatic vascular systems, J.Theoret. Biol. 355(2014), 194-207.
School of Mathematics, Renmin University of China, Beijing, 100872, P. R. China
E-mail address : [email protected] Department of Applied Computational Mathematics and Statistics, University of Notre Dame,Notre Dame, IN 46556, USA
E-mail address ::