Martine Audibert

Marlène Guillon

Jacky Mathonnat

Cerdi

Pascal Ferdi Wp

 Audibert

J >>

The challenge of improving efficiency of Soum Health Centers in Mongolia -What data tell us for Soum Health Centers in five provinces?

Keywords: Efficiency, Data Envelopment Analysis, Double bootstrapping, Primary care, Mongolia

come

Introduction

The global context

Since the mid-90's, Mongolia has made significant progresses in health. Between 1995 and 2014, life expectancy at birth increased from 61 to 69 years while under-five mortality rate and maternal mortality ratio were divided by four, from 85 to 24 per 1000 live births and from 205 to 46 per 100.000 live respectively (World Health Organization -WHO -estimates). They have been favored by various reforms, a radical change in the structure of economic growth and a favorable international environment.

But things changed, and new financial constraints have emerged. The economic situation of Mongolia has deteriorated considerably since 2012 up to 2016. At the root causes of the crisis, there are reasons coming from exogenous factors to which the authorities responded with measures that were not compatible with the precariousness of the situation. External demand weakened due to a continued dampening of the commodity market and slower growth in China, leading to a drop in Mongolian exports and fiscal revenues. There is therefore an extremely challenging environment for public finances while the state budget contributed 62% of the financing of health expenditure in 2016 (Center for Health Development, 2016).

The strong constraints linked to the macroeconomic situation also have repercussions on the financing of health by the insurance system. Health insurance is mainly financed by the contributions of employers and employees, and by the self-employed in the various sectors, including the informal sector of which herders constitute a large part. The insurance premiums subsidized by the state budget represent in 2016 nearly 16% of the total resources of the health insurance fund (Center for Health Development, 2016). And every year, because of economic difficulties a number of households, including herders and in-migrants, cannot pay their premiums and are therefore without health insurance (Ministry of Population Development and Social Protection, 2013).

As a corollary of this situation, very heavy constraints weigh on the public financing of health expenditures.

In addition to the fiscal constraint, the rise of non-communicable diseases and the sharp increase in private health expenditures are two major trends that underline the importance of improving efficiency in the health system. Non-communicable diseases, including cardiovascular diseases and diabetes, are the leading causes of years of life lost in Mongolia. Diabetes, which was quite rare in Mongolia 30 years ago, now hits about one in seven persons. The diabetes rate has been multiplied by six and diabetes-related mortality by four between 2005and 2014(WHO, 2015)). Regarding health expenditures, between 1995 and 2012, out-of-pocket (OOP) health expenditures rose in Mongolia from 12% to 42% of total health expenditures (WHO). The current level of OOP health expenditures is slightly above the average (39%) of the lower middle income countries (World Bank). Furthermore, global worldwide experience suggests that universal health coverage is difficult to achieve if OOP expenditures represent a high share of total health expenditures (Savedoff et al., 2012). Let's add that the Mongolian government has committed itself in 2013 to provide effective universal primary health care services for all and wants also to expand the costly Child Money Program to reach at least 80% of children from the poorest families .

In this very restrictive context for financing an ambitious health care policy, improving efficiency is a compelling challenge 1 .

A lack of knowledge about the efficiency of health systems in the countries of former

Soviet Central Asia.

The interest to study the efficiency of health facilities in Mongolia is reinforced by the fact that Central Asian post-soviet health systems have been ranked amongst the least studied over the world (McKee et al., 2012) despite the radical change in health care organization they experienced after the fall of the Soviet Union. As other former soviet health systems, the Mongolian health system was organized according to the Semashko model until 1990. The health system was then characterized by its centralized planning, administration and financing. It provided universally accessible and free basic health care but was heavily focused on curative and inpatient care while neglecting preventive and primary outpatient care. This health system became financially unsustainable after losing the subsidies of the Soviet Union and Mongolia engaged in major health 1 In 2017, the economic situation of Mongolia improved significantly due to favorable external environment and a sound macroeconomic policy (IMF, 2018). Mongolia benefited from a sharp increase in coal exports, positive terms of trade and large foreign investments in the commodities sector. This buoyant external context has been accompanied by a significant reduction in the fiscal deficit (increase in fiscal resources, containment of expenditures) and by a prudent monetary policy. But beyond this improvement, significant risks and challenges still lie ahead, including the level of debt, questioning regarding the dynamics of copper and coal exports and therefore the growth of fiscal revenues, limited diversification of the economy, and domestic political situation leading to ask about the political feasibility of the vigorous reform agenda in which the government is committed with a possible risk of skidding off due to popular opposition. In other words, that means that the current improvement in the macroeconomic situation does not detract from the strong need to improve efficiency in the health system. system reforms after 1990. These reforms focused on increasing efficiency through scaling back secondary and tertiary hospitals and strengthening primary care.

In recent years, a growing literature has emerged on the efficiency of primary health care facilities in in low-and middle-income countries. Available studies focus on African countries such as Ghana (Alhassan et al., 2015;Novignon and Novignon, 2017), Sierra Leone (Kirigia et al., 2011), Ethiopia (San Sebastian and Lemma, 2010) or Burkina Faso (Marschall and Flessa, 2009), on South-American countries such as Guatemala (Hernandez and Sebastian, 2014) or Chile (Ramírez-Valdivia et al., 2011) and on Asian countries such as Pakistan (Razzaq et al., 2013) or China (Audibert et al., 2013;Cheng et al., 2016). However, to the best of our knowledge, no published study has investigated with quantitative methods the efficiency of primary care facilities in countries with former Semashko health systems that are still facing major challenges linked to their over-reliance on inpatient and specialized care.

Soum Health Centers in the Mongolian health system

Mongolia is administratively divided into 21 aimags (provinces) which are in turn split, for rural areas, into 329 soums (districts). The Mongolian health system is based on a two-tier model that provides primary health care at the primary level and specialized health care at secondary and tertiary levels (Tsilaajav et al., 2013).

In rural areas, primary health care is mainly delivered by Soum Health Centers (SHCs) on which our analysis focuses. SHCs are public health facilities owned by local governments and run by salaried staff, i.e. doctors, nurses, midwives and support staff. SHCs provide antenatal and postnatal care, minor surgeries and normal deliveries as well as preventive activities such as immunization, diabetes/hypertension testing or health education. In 2011, a revision of the health act stressed the focus of SHCs on primary health care and prevention while reducing their hospital functions.

Nevertheless, SHCs kept providing both outpatient and inpatient services. The double mission of SHCs is explained by the geographic and demographic characteristics of Mongolia where the rural and nomadic populations are sparsely distributed. In SHCs, inpatient admissions, outpatient visits or routine immunization are provided free of charge for all patients. However, the reimbursement rate of outpatient essential drugs varies between 5% and 91% of cost price. Only patients with Social Health Insurance (SHI) entitlement, i.e. patients who paid their insurance contributions for the last 12 months, are eligible for outpatient drug reimbursement. An important source of OOP expenditures for SHCs' patients is then copayment for essential drugs or full drug costs if uninsured. In 2011, pharmaceuticals represented 94% of health OOP payments among the very poor in Mongolia (Tsolmongerel et al., 2011).

The paper is structured as follows. Section 2 details the methodology and Section 3 presents the dataset. Section 4 displays the results that are discussed in Section 5.

Methods

We use a two-stage model to estimate SHCs' efficiency and its determinants. We apply the double bootstrap DEA procedure developed by Simar and Wilson (2007) to calculate bias-adjusted DEA scores and to study the factors associated with SHCs' efficiency.

Technical efficiency refers to the capacity of a Decision Making Unit (DMU), a SHC in our case, to transform a certain amount of inputs into a certain quantity of outputs through the production process. In an output orientation, technical efficiency reflects the ability of a SHC to obtain a maximal level of outputs from a given set of inputs. In an input orientation, it reflects the ability of a SHC to minimize inputs given a level of production. For the efficiency analysis, as SHCs have to deal with fixed amounts of resources, the output orientation was chosen. Two methods are available to estimate the technical efficiency of SHCs: a parametric one using the Stochastic Frontier Analysis (Aigner et al., 1977) and a non-parametric one, which includes several methods such as the DEA (Farrell, 1957;Charnes et al., 1978) or the Free Disposal Hull (FDH) (Deprins et al., 1984). As the parametric method requires an ad hoc assumption on the functional form of the production function we use a non-parametric approach. Among the different non-parametric methods, DEA is the most often used to estimate health facilities' efficiency (Hollingsworth, 2008). The purpose of the DEA method is to construct a piecewise linear envelopment frontier over the data points such that all observed points lie on or below the production frontier. The efficiency level of each SHC is then measured by calculating the difference between the point representing the observed values of inputs and outputs of this SHC and the frontier that determines the best observed practices in the whole sample. In DEA, technical efficiency can be measured using either constant or variable returns to scale. The assumption of constant returns to scale is appropriate if all SHCs operate at the optimum scale. If it is not the case, efficiency estimates will include scale economies while we aim to measure pure technical efficiency. We therefore consider variable returns to scale in the DEA analysis.

The use of two-stage models to study the factors associated with efficiency has been challenged (Simar and Wilson, 2007). The controversy lies on the nature of the efficiency scores used as the dependent variable in the second-stage. The instrumentalist approach considers the estimated efficiency scores as descriptive measures of the relative technical efficiency of sampled DMUs (McDonald, 2009;Ramalho et al., 2010). Efficiency scores can then be treated as a standard dependent variable and classical regression analyses carried out in the second stage are thought to provide valid inference. On the contrary, the conventionalist approach looks at efficiency scores as estimates of true efficiency scores (Simar and Wilson, 2007;Simar and Wilson, 2013). Efficiency scores are assumed to reflect not only efficiency but also noise because of measurement errors in inputs and outputs. Moreover, since the estimated frontier can only be a subset of the true but unknown technology, efficiency scores are biased due to sampling variability. Additionally, given that efficiency scores estimated through non-parametric approaches are dependent on each other, efficiency estimates are considered serially correlated. Consequently, the conventionalist approach regards the standard approaches of inference used in the second-stage regression analysis as invalid. Simar and Wilson (2007) propose a double bootstrap procedure that produces bias-adjusted DEA score and provides valid inference for the regression of non-parametric DEA efficiency estimates on a set of explanatory variables in the second stage. First, bias-corrected DEA scores are calculated by subtracting the bootstrap bias estimate from the original DEA estimates. Bias-corrected DEA scores are then regressed on a set of explanatory variables using a bootstrapped truncated regression that produces unbiased parameter estimates and confidence intervals. We use the Algorithm 2 of Simar and Wilson (2007) after conducting outliers' detection based the method developed by Simar (2003). Algorithm 2 proceeds in the following seven steps: S1: Estimate output-orientated Shephard DEA efficiency scores with varying returns to scale for all SHCs in the sample.

For each SHC i=1,…,n with observed outputs and inputs , output-oriented Shephard DEA efficiency estimate is the solution of the following linear program:

=(max{θ > 0 | θ ≤ ∑ ; ≥ ∑ ; ∑ = 1 ; ≥ 0, i = 1, … , n}) (1)
For all i=1,…,n, ≤ 1. A SHC is technically efficient if =1 and technically inefficient if < 1.

The proportional increase in output that can be achieved by an inefficient SHC with no increase in inputs is equal -1 *100 per cent.

S2:

Estimate Equation 2 by employing truncated maximum likelihood to yield estimates β and .

= β + ≤ 1 (2)
In Equation (2), is a vector of variables assumed to influence the values of and , β is a vector of parameters to be estimated and is a continuous iid random variable independent of with distribution N(0, ).

S3:

For each i = 1,…,n, loop over the following four steps (i-iv) times (100 times in our case) to obtain a set of bootstrap estimates = { * } .

(i) For each i=1,…,n, draw from the N(0, ) distribution.

(ii) For each i=1,…,n, compute * = + .

(iii) Construct a pseudo data set ( * , * ) where * = and * = ( * )

Ferdi WP n°226  Audibert, M., Guillon, M., and Mathonnat, J. >> The challenge of improving efficiency … 6 (iv) Using the pseudo data set and Equation 1, compute pseudo efficiency estimates * for all i=1,…,n.

S4:

For each SHC i=1,…,n, compute the bias-corrected estimator = -Bias( ) where the bias term is estimated as: Bias( )=( ∑ * )-.

S5:

Use the method of maximum likelihood to estimate the truncated regression of on to yield estimates and .

S6: Loop over the following three steps (i-iii) times (1000 times in our case) to yield a set of bootstrap estimates Φ = {( * , * ) } .

(i) For each i=1,…,n, draw from the N(0, ) distribution.

(ii) For each SHC i=1,…,n, compute * * = + .

(iii) Use the maximum likelihood method to estimate the truncated regression of * * on to yield estimates * and * .

S7:

Use the boostrap values in Φ and the original estimates and to construct estimated confidence intervals for each element of β and . The (1-α) per cent confidence interval of the j th element of vector β, where 0<α<1, is constructed as the Pr(-

/ ≤ * -≤ -/ ) ≈ 1- α such that the estimated coefficient is [ + / * , + / * ].

Data

The analysis focuses on five regions (aimag) for which we obtained complete and reliable data on SHCs' activity and resources use between 2013 and 2015: Arkhangai, Bulgan, Khuvsgul, Sukhbaatar and Zavkhan (Figure 1). nurses might not substitute to doctors in the realization of all medical acts. In the main analysis (model 1), we also consider three outputs: the number of preventive visits, the number of nonpreventive visits and the number of inpatients. We differentiate preventive and non-preventive visits since the 2011 revision of the Mongolian health act stressed out the focus of SHCs on primary preventive care. The number of inpatients is included since SHCs also have for mission to provide inpatient care for rural populations whose access to secondary level health facilities (aimag general hospitals) in aimag centers is difficult. As a robustness analysis, we consider alternative outputs in models 2 to 8 (Table 1). We successively decompose non-preventive outpatient visits into SHC and home visits, replace the number inpatients by the number of beds days and substitute the number of preventive visits by the number of controlled pregnant women and the number of diabetes and hypertension tests performed. Non-preventive visits (home)

Potential determinants of efficiency

We consider two kinds of variables that might affect SHCs' efficiency; external variables that are not under the control of the SHCs; and internal variables whose value is discretionary to the health authorities.

Regarding external variables, we first consider the characteristics of the soum population: the size of the population, the poverty rate measured by the poverty headcount ratio (Coulombe and Altankhuyag, 2012) 2 , the age structure of the population measured by the dependency ratio and the share of the population living in ger. We expect to find a positive correlation between the population of the soum and the efficiency of SHCs through a demand effect. The poverty rate in the catchment area of a SHC might have several effects on its efficiency. If poorer people tend to exhibit more deteriorated health states, a higher poverty rate in the catchment area of a SHC could lead to an increase in its activity and efficiency. On the other hand, given the existence of OOP payment for outpatient drugs in SHCs, a higher poverty rate might increase health care renouncements. Furthermore, travelling expenses to the SHC might constitute a significant financial burden for poor population, leading to more health care renouncements. The age structure of the population in its catchment area may also affect the activity of a SHC. As children and seniors consume more of both preventive and curative care (Gan-Yadam et al., 2013), a high dependency ratio in the soum might be associated with more activity and increased efficiency for the SHC. The share of the soum population living in ger is also likely to influence SHCs' activity. As nomads live in remote areas, they might have difficulties in accessing health facilities. On the other hand, their harsher living conditions might lead to a higher prevalence of (infectious) diseases, which could translate into more demand for SHCs.

We also consider two geographical variables, the area of the soum in square kilometers and the distance of the SHC to the aimag general hospital in kilometers. The area of the soum is likely to affect the volume of demand faced by SHCs. A negative correlation between the area of the soum and the efficiency of SHCs could be observed if populations living in larger soums experience more difficulties to reach the SHC. The distance of the SHC to the aimag general hospital might be positively correlated with its activity and efficiency if the lack of accessibility of the soum population to the aimag general hospital translates into more demand for the SHC. Moreover, people might be less keen to bypass the referral system if the aimag general hospital is located further.

Based on the data available for the purpose of this study, we include one internal variable that might influence SHCs' efficiency, the share of doctors in medical staff (in %). The share of doctors might be positively correlated with SHCs' efficiency if doctors are more efficient than nurses in the delivery of health care and then treat more patients. However, this could also lead to fewer revisits through a "quality of care" effect, and then to a decrease in the global level of SHCs' activity and efficiency.

Main characteristics of SHCs and of their catchment area

Table 2 provides key statistics for the sample of SHCs used in model 1.

The mean number of inpatients is 331 while the mean number of bed days is 2361. In average, SHCs realized 3119 preventive outpatient visits and 4593 non-preventive outpatient visits per year.

Non-preventive visits in the SHC are higher than home visits with means of 3221 and 1373 respectively. Among SHCs of our sample, the mean annual number of controlled pregnant women is 401 while the mean annual number of hypertension and diabetes tests is 843. In average, SHCs have 9.7 beds. The mean numbers of doctors and nurses are 2.6 and 7.2 respectively. In average, 39.2% of the population lives below the poverty line and 77.6% lives in ger. The correlation between the poverty headcount and the share of the population living in ger is low (0.0409) and non-significant at a 10% level. The mean dependency ratio is high and equal to 51.8%. Soums measure an average of 4061 square kilometers and host an average of 3199 habitants. The mean distance of SHCs to the aimag general hospital is 144 kilometers and the mean share of doctors in medical staff is 27.1% among SHCs of our sample.

Results

As DEA scores are highly sensitive to the presence of outliers in the sample, we run a preliminary analysis to detect outliers. We use the leave-one out outliers' detection methodology developed by Simar (2003). Details of the procedure are available in Appendix A. After several repetitions to control for masking effects, we exclude three outliers: Tynel and Byrentogtokh SHCs in Khuvsgul in 2013 and Ikh-Uul SHC in Zavkhan in 2013.

Efficiency

Factors having an impact on efficiency

Table 5 presents the results of the truncated regressions using bias-adjusted DEA scores. We first present the results for the main model (model 1) and then discuss the heterogeneity of results by models.

Looking at model 1, we find a negative and significant correlation between the poverty headcount and the efficiency score. This result is confirmed by the use of two alternative measures of poverty;

the poverty gap index and the poverty severity index (results not shown). Results of the regression analysis for model 1 show a positive and significant association between the dependency ratio and the efficiency of SHCs. Then, the presence of more children and elders in its catchment population enhances the efficiency of a SHC. The size of the soum population is also positively correlated with the efficiency score. This indicates that the low efficiency exhibited by some SHCs is related to the small demand for medical care they face in sparsely populated areas. On the contrary, we find no correlation between the area of the soum and the efficiency of SHCs once controlled for the population size. The share of the soum population living in ger is positively correlated with the efficiency of SHCs. This implies that SHCs facing more nomadic populations exhibit better efficiency. The proximity of SHCs to the aimag general hospital does not impact their efficiency.

Moreover, we find no significant correlation between the share of doctors in medical staff and the efficiency of SHCs.

The sign and significance of correlations are unchanged for several variables whatever the model considered. This is the case for the poverty headcount, the share of the catchment population living in ger, the distance to the intersoum hospital and the size of the population (except in model 4). On the contrary, the dependency ratio loses statistical significance when preventive visits are replaced by the number of controlled pregnant women and the number of diabetes/hypertension tests. The area of the soum is positively correlated with SHCs' efficiency only in models where bed days are used instead of inpatients. Finally, the share of doctors in medical staff is negatively correlated with SHCs' efficiency in models 3, 5 and 7 while it is not significant in other models.

Discussion

We study the efficiency of SHCs which provide primary health care for rural populations in Mongolia. Between 2013 and 2015, SHCs of our sample exhibit a low mean efficiency (0.682)

indicating that health care production could, in average, be increased by 46.67% without an increase in inputs. Our results point to a decrease in efficiency of SHCs over time, in particular between 2014 and 2015, and to a heterogeneity of SHCs' efficiency across regions.

Our results show that demand-side factors impact SHCs' efficiency. First, the size of the soum population affects positively the efficiency of SHCs because SHCs located in more populated areas face a higher demand that spurs their activity. This result is in line with those found for rural primary health centers in Chile (Ramírez-Valdivia et al., 2011) or Burkina Faso (Marschall and Flessa, 2009), where more population in the catchment area is associated with better efficiency. However, this result diverges from those of San Sebastian and Lemma ( 2010 (2013), who find a negative correlation between rural income and township hospitals' efficiency in China. Then, providing free outpatient drugs in SHCs, even for those not covered by the SHI, could increase SHCs' activity and efficiency. Besides ensuring free primary health care, access to care of poor populations and SHCs' efficiency could benefit from the creation of better and low cost transportation infrastructures in the long-run or -which would be more feasible considering the huge surface of Mongolia -to benefit from transport costs partially subsidized by the SHI.

SHCs' efficiency is positively associated with the share of the soum population living in ger. Despite a more difficult access to health facilities, people living in ger experience harsher living conditions that may deteriorate their health states and lead them to demand more medical care than sedentary inhabitants. This could explain why SHCs facing more nomadic inhabitants achieve better efficiency. The age structure of the catchment population also impacts SHCs' efficiency. We find that SHCs facing more children and seniors exhibit higher efficiency, likely because these population categories consume more preventive and curative health care than middle-age adults (Gan-Yadam et al., 2013). However, the opposite result is found by Ramírez-Valdivia et al. (2011) in

Chile where the shares of the population under 6 and over 65 are negatively associated with the efficiency of primary health care practices in rural municipalities.

Our results also point to the lack of association between SHCs' efficiency and the distance to the aimag general hospital. Thus, the activity and efficiency of SHCs are not affected by the proximity of a higher referral level heath facility. This indicates that the self-referral of patients to the aimag general hospitals is not a major issue for the rural health system in Mongolia. Two previous studies, conducted in Ethiopia (San Sebastian and Lemma, 2010) and China (Audibert et al., 2013), also find no association between primary health facilities' efficiency and their distance to the closest higher level health center.

SHCs' efficiency is positively correlated with the area of the soum only when considering the number of bed days instead of the number of inpatients. To understand this result, we ran additional regression analyses3 . Controlling for soum attributes and population characteristics, we find that the soum area is positively and significantly correlated with both the number of bed days and the mean length of stay. However, this variable is not associated with the number of inpatients. This result might indicate that patients living further from the SHC benefit from longer inpatient stay.

With the exception of model 1, the share of doctors in medical staff is negatively correlated with SHCs' efficiency when considering inpatients rather than bed days. Further regression analyses4

show that the share of doctors is positively correlated with the number of bed days and the mean length of stay, but not associated with the number of inpatients. Then, the negative association between the share of doctors and efficiency, when considering the number of inpatients as an output, might be linked to longer inpatient stays decided by doctors compared to nurses. These longer inpatient stays might be associated with higher quality of care and fewer revisits.

Our study is not without limitations. Data have been collected at the SHC level and we could only gather complete and reliable data for SHCs of 5 out of 20 Mongolian rural regions. Then, our results regarding the determinants of SHCs' efficiency cannot be generalized to the entire country as the partial data (or lack of data) for the SHCs of other regions did not make it possible to integrate them into the analysis. Moreover, the analyses made in this study to highlight the factors that contribute to explain the differences in SHCs' efficiency were conditional on available data relevant for the purpose of this study. That is, there are several potential explanatory fields -which can be important from the point of view of the regulation of the health system -that could not be taken into consideration. This is the case of the role of financial incentives which are complementary to remunerations and whose effect on efficiency, positive or negative, is a priori indeterminate.

Further studies should integrate such variables in order to provide more guidance to the Mongolian health authorities from a policy-oriented perspective. Despite these limitations, our study can provide useful information for policy makers and constitutes the first attempt to measure health care facilities' efficiency and efficiency determinants in central Asian post-soviet health systems. Our results are likely to be of interest for other countries that transitioned from the Semashko health system after the fall of the Soviet Union but that are still heavily relying on inpatient and specialized care, especially those sharing geographic characteristics with Mongolia such as the low population density in rural areas.

Figure 1 :

 1 Figure 1: Regions of provenance for SHCs data

  ) orAudibert et al. (2013). The former find no association between the sub district's population and the efficiency of rural health posts in Ethiopia while the latter find no correlation between the population density and the efficiency of township hospitals in China. A logical policy implication of this result would be -taking careful consideration of the characteristics of the local context -to reduce the resources of SHCs which face a too low demand. Indeed, our analyses indicate that most SHCs could achieve similar activity levels with reduced resources, therefore not affecting equity of access to health care between rural and urban populations.

Facing

  poorer populations tend to lessen SHCs' efficiency, most likely because OOP health expenditures or transportation costs induce health care renouncements among low-income populations. This result supports the findings ofDorjdagva et al. (2017), who point to pro-rich inequalities in the use of inpatient services in SHCs. This is also in line with results of Audibert et al.

Table 1 : Presentation of the alternative models

 1

	Outputs

Table 2 : Descriptive statistics of the sample

 2

		Obs	Mean	Std. Dev.	Min	Max
	Inpatients	221	330.6968	108.5792	108	832
	Preventive visits	221	3119.018	2147.682	116	13980
	Non-preventive visits	221	4593.294	2571.983	775	15755
	Bed days	220	2360.736	742.8844	971	5317
	Non-preventive visits in SHC	221	3221.176	2007.224	553	13118
	Home non-preventive visits	221	1372.692	932.167	114	4979
	Pregnant women controlled	221	401.3484	183.4399	104	1110
	Hypertension and diabetes tests	204	843.2255	493.0895	0	2458
	Beds	221	9.642534	3.64365	4	21
	Nurses	221	7.162896	2.748666	3	15
	Doctors	221	2.556561	1.141257	1	7
	Poverty headcount	221	39.2009	6.645155	27.4	50.5
	Dependency ratio	221	51.79909	4.620586	40.9	66.1
	Living in ger (%)	221	77.63484	16.13893	11.2	95.5
	Population	221	3198.91	1291.659	938	6549
	Area (km2)	221	4060.955	2592.972	626.9087	17002.7
	Distance general hospital	221	144.4344	74.64943	22	307
	Doctors (% medical staff)	221	27.13697	9.713719	6.25	60

Table 3

 3 to 0.62 in 2015 (t = 5.2232, p < 0.001). On the contrary, efficiency scores' heterogeneity does not evolve over time as the standard deviation is equal to 0.11 in all years. The mean efficiency score is between 0.665 and 0.748 in models 2 to 8. The highest correlation of efficiency scores between models is seen for model 1 and model 3 where the number of non-preventive visits is split between home and SHC visits. The correlations between the efficiency scores of model 1 and alternative models are much lower when the number of preventive visits is replaced by the number of controlled pregnant women and the total number of diabetes/hypertension tests. This result is logical since preventive visits in SHCs integrate other activities such as immunization or health education.

gives the summary statistics of the efficiency scores for all models. Mean efficiency in the main model (model 1) is equal to 0.682, which means that, in average, SHCs could increase output by 46.67% without an increase in inputs. Mean efficiency decreases over time in model 1, from 0.72 in 2013

Table 3 : Summary statistics of the efficiency scores

 3 Table 4 displays the efficiency scores by region and by year for model 1. The mean efficiency score decreases in all regions between 2013 and 2015, though at different rates. The decrease in mean efficiency is the highest in Khuvsgul with a drop of -19.71% between 2013 and 2015. Thereby, Khuvsgul ranks as the most efficient region in 2013 but as the second least efficient in 2015. The lowest decrease in mean efficiency is observed in Bulgan (-5.08%). Consequently, Bulgan is the second most efficient region in 2015 while being the least efficient in 2013. Between 2013 and 2014, mean efficiency increases in Bulgan and Sukhbaatar, decreases in Khuvsgul and Zavkhan and remains stable in Arkhangai.

		Observation	Mean	Standard deviation	Min	Max	Correlation with model 1
	Model 1	221	0.6817818 0.1204662	0.3889395	0.9368711	-
	Model 2	220	0.6652675 0.1322318	0.3462393	0.9307621	0.8747***
	Model 3	221	0.7024325 0.1210076	0.3994593	0.9355465	0.9502***
	Model 4	220	0.686424 0.1345401	0.3630884	0.9390115	0.8187***
	Model 5	200	0.7270786	0.12205	0.366824	0.9329448	0.6637***
	Model 6	199	0.7069749 0.1297229	0.333037	0.935676	0.5606***
	Model 7	200	0.7471827 0.1208079	0.3751675	0.945289	0.6686***
	Model 8	199	0.7212512 0.1289889	0.3379224	0.9347335	0.5626***

Table 4 : Efficiency scores by region and year 2013-2015 2013 2014 2015 2013-2015 variation Mean score Ranking Mean score Ranking Mean score Ranking Mean score Ranking

 4

	Arkhangai	0.6662246	4	0.6947109	3	0.6934461	5	0.6105169	3	-12.12%
	Bulgan	0.6737065	3	0.6695225	5	0.7157627	3	0.6355126	2	-5.08%
	Khuvsgul	0.6870264	2	0.75905	1	0.7095282	4	0.6094476	4	-19.71%
	Sukhbaatar	0.6588835	5	0.6845962	4	0.717495	2	0.5745593	5	-16.07%
	Zavkhan	0.7101955	1	0.7524139	2	0.7385126	1	0.6441042	1	-14.39%

Table 5 : Results of truncated regression analyses Bias-adjusted DEA score

 5

		Model 1	Model 2	Model 3	Model 4	Model 5	Model 6	Model 7	Model 8
	Poverty headcount	-0.00327**	-0.00531***	-0.00430***	-0.00538***	-0.00336***	-0.00286**	-0.00445***	-0.00369**
		(0.00132)	(0.00142)	(0.00132)	(0.00141)	(0.00125)	(0.00146)	(0.00138)	(0.00145)
	Dependency ratio	0.00808***	0.00470**	0.00768***	0.00456**	0.000144	-0.00375*	0.000467	-0.00338
		(0.00197)	(0.00206)	(0.00196)	(0.00210)	(0.00201)	(0.00216)	(0.00204)	(0.00228)
	Living in ger (%)	0.00169***	0.00218***	0.00197***	0.00247***	0.00127***	0.00141***	0.00160***	0.00184***
		(0.000468)	(0.000503)	(0.000462)	(0.000514)	(0.000460)	(0.000520)	(0.000490)	(0.000534)
	Population	0.0000152**	0.0000132*	0.0000166***	0.0000109	0.0000662*** 0.0000624*** 0.0000635*** 0.0000628***
		(0.00000619) (0.00000695)	(0.00000614) (0.00000678) (0.00000677) (0.00000746) (0.00000667) (0.00000773)
	Area (km2)	-0.00000203 0.00000910*** -0.00000154 0.0000115*** 0.000000899 0.00000695** 0.000000246 0.00000721**
		(0.00000283) (0.00000327)	(0.00000281) (0.00000347) (0.00000292) (0.00000353) (0.00000313) (0.00000356)
	Distance general hospital -0.0000194	0.00000962	-0.0000517	-0.0000138	-0.0000768	-0.0000798	-0.0000414	0.0000176
		(0.000106)	(0.000116)	(0.000105)	(0.000117)	(0.000113)	(0.000129)	(0.000114)	(0.000126)
	Doctors (% medical staff)	-0.00125	-0.0000812	-0.00160*	-0.000503	-0.00289***	-0.00178	-0.00281***	-0.00164
		(0.000898)	(0.000967)	(0.000890)	(0.00100)	(0.000934)	(0.00109)	(0.000965)	(0.00106)
	Year = 2014	0.0400**	0.00284	0.0417**	0.00418	0.0258	-0.0141	0.0221	-0.0258
		(0.0198)	(0.0217)	(0.0197)	(0.0213)	(0.0196)	(0.0222)	(0.0207)	(0.0222)
	Year = 2015	-0.0702***	-0.115***	-0.0741***	-0.123***	-0.0327*	-0.0720***	-0.0385**	-0.0804***
		(0.0181)	(0.0203)	(0.0180)	(0.0206)	(0.0183)	(0.0196)	(0.0194)	(0.0208)
	Constant	0.269***	0.423***	0.337***	0.448***	0.626***	0.764***	0.657***	0.747***
		(0.101)	(0.107)	(0.0999)	(0.112)	(0.0985)	(0.110)	(0.100)	(0.111)
	N	221	220	221	220	200	199	200	199
	Standard errors in parentheses ; * p<0.10, ** p<0.05, *** p<0.01

Poverty data are not available yearly at the soum level in Mongolia. Therefore, we use the United Nations Development Programme data for 2011, which are the latest poverty data available at the soum level.

Available from the authors upon request.

Idem.

* This study is part of a joint research project FERDI and School of Public Health, Mongolian National University of Medical Sciences (MNUMS), under the coordination of Jacky Mathonnat. It has been supported by the Agence Nationale de la Recherche of the French government through the program ''Investissements d'avenir ANR-10-LABX-14-01'' , by the Foundation for Studies and Research on International Development (FERDI) and by the School of Public Health, MNUMS.

The funding sources had no role in the study design, in the collection, analysis and interpretation of data nor in the writing of this report. The authors thanks their colleagues L. Munkh-Erdene, L. Gantugs and A. Purevgered (School of Public Health) for their fruitful collaboration, as well as the medical teams of the sum health centers and the administrative officials met during the field missions. They also thank Dr Laurène Petitfour (CERDI) for her involvement in a preliminary phase of the project.

• Aigner, D., Lovell, C. K., and Schmidt, P. (1977). Formulation and estimation of stochastic frontier production function models. Journal of econometrics, 6(1), 21-37. • Alhassan, R. K., Nketiah-Amponsah, E., Akazili, J., Spieker, N., Arhinful, D. K., and de Wit, T. F. R. (2015). Efficiency of private and public primary health facilities accredited by the National Health Insurance Authority in Ghana. Cost Effectiveness and Resource Allocation, 13(1), 23.

• Audibert, M., Mathonnat, J., Pelissier, A., Huang, X. X., and Ma, A. (2013)


Laurent Sagart

Crlao

OC *-r in early Chinese loans to Bùyāng, and related issues

 reconstruct a single liquid coda in Old Chinese, noted *-r. This *-r evolves to MC *-n or *-j depending on dialect. In certain Tai and Mon-Khmer languages, some Chinese loanwords with the OC *-r coda end in -i, even when the MC pronunciation has -n. We pointed out that "If the original coda had been [n], there is no reason it could not have been borrowed with [n]". 1 The fate of OC *-r in Chinese loans to Bùyāng 布央, and in Càijiā 蔡家, is similar.

Layers of Chinese loanwords in Lángjià Bùyāng

Bùyāng 布央 is a Kra-Dai language spoken in Yúnnán. It belongs to the Kra (aka 圪央 Ge-yang) branch. Its Lángjià 朗架 dialect was described by [START_REF] Jǐnfāng | Buyang Yu Yanjiu 布央語研究[END_REF]. Bùyāng has borrowed extensively from Chinese since first contact between Chinese and Kra-Dai, which presumably occurred when Qín Shǐ Huáng's armies established a Chinese state in the region, around 200 BCE. A preliminary analysis of the stratification of Chinese loanwords based on the reflection of Chinese tones and segments points to at least six layers of direct or indirect Chinese loanwords.

Treatment of Old Chinese *-r in loans to Bùyāng

'dog'

One might be tempted to dismiss the comparison between 犬 quǎn 'dog' and Bùyāng qɔi 24 as phonetically adventurous. Yet absence of aspiration in qɔi 24 is expected in an early Chinese loanword, as no aspirated stops can be reconstructed at the earliest layer of KD; the uvular initial is regular for a Chinese type-A word borrowed before Late Middle Chinese, whether the OC initial vas velar or uvular; and the tone matches Chinese shǎngshēng in a Layer-II word.

'dog ' in Càijiā.

We have collateral from Càijiā 蔡家, an archaic Chinese dialect of Guìzhōu, for an *-r coda in 犬 'dog'. Càijiā was recently investigated by Lü Shanshan and Miguel Cortiço dos Santos. They were generous to share with me the results of their fieldwork. They recorded 'dog' as khwi 55 . Tone 55 is the main Càijiā reflex of OC -ʔ, suggesting we are dealing with the Chinese word, rather than with a loan from a neighboring ST language. With one counterexample in the Càijiā 200-word-list at my disposal (ȵi³³ 'person', OC *niŋ), Càijiā reflects OC *-in and *-iŋ, perhaps also OC *-en, as -ɪɴ (/iN/) or -iɴ (/jiN/):

ȵiɴ³¹ 'year' 年 *C.nˤi[ŋ] kʰɪɴ 55 'day' 天 *l^ˤi[n] ɕɪɴ³³ 'new' 新 *C.si[n] mɪɴ³¹ 'people' 民 *mi[ŋ] kʰɪɴ³³ 'to pull' 牵 *[k]ʰˤi[n] kɪɴ³³ 'to see' 見 *[k]ˤen-s

Correspondence of OC *-r in Written Tibetan

. [START_REF] Hill | Cognates of Old Chinese *-n, *-r, and *-j in Tibetan and Burmese[END_REF] What happens to ST words with high vowels followed by the coda *-r in other ST languages ?

In Table 8

Conclusion: a third candidate for a non-Chinese phonological innovation

Thus it appears that after high vowels, the PST coda *-r was maintained in Chinese but changed to *-j, merging with original *-j in all Sino-Tibetan languages outside of Chinese. This constitutes another candidate for a shared innovation of non-Chinese ST. Another candidate was proposed in [START_REF] Handel | Rethinking the medials of Old Chinese: Where are the r's?[END_REF], yet another in [START_REF] Sagart | A candidate for a Tibeto-Burman innovation[END_REF].

  examines the correspondences in Tibetan of OC *-r as reconstructed by[START_REF] Baxter | Old Chinese: a new reconstruction[END_REF]. He finds two: one to Tibetan *-r, for which he reconstructs *-r in the parent language (I use the term 'Proto-Sino-Tibetan'); and another to Tibetan *-l, for which he uses the index *-rl. If we now assign *-r to the Old Chinese word for 'dog' we seem to need a third correspondence: OC *-r to Tib. -zero, to account for Tib. ཁ khyi 'dog'. Three comparisons are known to me which exhibit that correspondence: dog 犬 *[k]ʷʰˤ[e]rʔ > khwenX > quǎn 'dog' Tib. ཁ khyi water 水 *s.turʔ > sywijX > shuǐ 'water; river' Tib. ཆ chu lip 脣 *sə.dur > zywin > chún 'lip' Tib. མཆ mchu Table 6. OC *-r to Tibetan zero after an OC high vowelHill's evidence for the correspondence OC *-r = Tib. *-r involves only words with nonhigh vowels, *a, *e, *ə, *o on the Chinese side (the Old Chinese reconstructions and Tibetan transcriptions have been updated; the evidence is limited to words B&S reconstruct with *-r or *-[r]): distribution between two Chinese-Tibetan sound correspondences 7. The correspondence OC *-r : Tib. zero beyond Chinese and Tibetan.

Table 1

 1

		Yīnpíng	Yángpíng	Yīnshǎng	Yángshǎng	Yīnqù	Yángqù
		陰平	陽平	陰上	陽上	陰去	陽去
	Early I		11	24	24 (some 11)		11
	Early II	54	312		24		11
	Zhuàng	24	11	11	11	54/31 54/31
	Late Middle Chinese	54	11	(24 ?)	?/31(QZ)	11	312
	Mandarin I	33	31	54	54/24(QZ)		24
	Mandarin II	11	24				11
	It is not entirely clear whether the Early I and II layers are chronologically distinct or whether they
	represent distinct contemporary sources. The Chinese twelve earthly branches, used as year names
	in Chinese imperial times, are part of the Early I layer:

. Layers of Chinese loanwords in Lángjià Bùyāng (in the Yángshǎng column, a contour followed by the mention '(QZ)' is conditioned by a voiced obstruent initial).

Table 2 :

 2 Borrowed year names in BùyāngBùyāng has also borrowed extensively from Zhuàng 壯, the locally dominant Kra-Dai language, and a member of the Tai subgroup. As Zhuàng loans themselves include much Chinese material, a Zhuàng layer is included in Table1above. Chinese loans to Zhuàng were made over as long a

period as have loans to Bùyāng: material borrowed by Bùyāng from Zhuàng was potentially borrowed at the time of the earliest contact with Chinese. The Late Middle Chinese and Mandarin layers were borrowed more recently.

Table 3 :

 3 Bùyāng -i for Old Chinese *-r > MC -n

	The first three words in Table 3 are reconstructed with OC *-r, and the fourth, with *-[n], standing
	for either with *-n, *-ŋ or *-r, with *-n the preferred option.

Four words from the Early and Zhuàng layers which end in -n in Middle Chinese have -i in Bùyāng: BS 2014 Bùyāng old Chinese loans layer 辰 *[d]ər > dzyin > chén '5th earthly branch' ɕi 11 'year of the dragon' Early I 寅 *[ɢ](r)ər > yin > yín '3rd earthly branch' ŋi 11 ' year of the tiger' Early I 蒜 *[s]ˤor-s > swanH > suàn 'garlic' ɕu:i 54 'garlic' Zhuàng 犬 *[k]ʷʰˤ[e][n]ʔ > khwenX > quǎn 'dog' qɔi 24 'dog' Early II

The normal Bùyāng treatment of OC *-n appears to be -n in the same layers.

Table 4

 4 An exception is 斤 *[k]ər > kj+n > jīn 'axe; catty', treated by Bùyāng as qaːn54 'ax', apparently in the Early II layer, despite having been reconstructed with *-r. This word cannot be part of the Late Middle Chinese layer, on account of its uvular initial. With this unexplained exception, the Bùyāng loan data are consistent with the view that the earliest Chinese donor(s) to that language distinguished between OC

*-r and *-n. This argues that the coda in the Old Chinese word for 'dog': 犬 *[k]ʷʰˤ[e][n]ʔ (BS 2014) was *-r.

Table 5 :

 5 Caijia -ɪɴ/-iɴ reflecting OC *-in and *-enIn general, with exceptions, 3 -ɴ is the expected Càijiā treatment of Old Chinese nasal endings *-m, *n, *-ŋ.

The treatment of OC *-r is different. Outside of a single example of a word reconstructed with OC *-r and ending in Càijiā -ɴ: 近 *N-kərʔ > gj+nX > jìn 'near', Càijiā tɕɪɴ³³, Càijiā loses *-r, even when its MC reflex is -n: 'dry' 乾, Càijiā ka 33 , OC *[k]ˤar; probably also in 'liver' 肝, Càijiā ka 33 , OC *s.kˤa

[r]

, and in 'to live' 存, Càijiā tə 31 , OC *[dz]ˤə

[n]

. Loss of OC *-r appears to be the normal Càijiā development. This implies, again, that the word for 'dog': khwi

55

developed out of an Old Chinese prototype with *-r coda.

  we add forms from Bodo, Lushai and Karen, representing the Sal, Kuki-Chin-Naga and Karen groups of ST, respectively.

	Tib.	OC (BS, revised)	Bodo (Bhat)	Lushai (Lorrain) Proto-Karen
					(LuangThongkum)
	ཁ khyi 'dog' 犬 *[k]ʷʰˤ[e]rʔ 'dog' səy(má) 'dog'	ui 'dog'	thwi B 'dog'
				tui 'water'	thej A 'water'
	མཆ mchu 'lip' 脣 *sə.dur 'lip'	(gusu)təy 'lip'	--	n.a.

ཆ chu 'water' 水 *s.turʔ 'water' dəy 'water'

Table 8 :

 8 An informal survey of forms for the etyma in Table8in other ST languages shows mostly j-type endings. Consequently as I argued here, the correspondence OC *-r : Tib. zero : other ST -i/-j probably reflects PST *-r after high vowels: 'water' and 'lip' from PST *-ur, 'dog' from PST *-ir. We would like to disambiguate the Old Chinese vowel to *[k]ʷʰˤirʔ, but under the bottom-up reconstruction principles adopted in Baxter and Sagart 2014, we must remain with *[k]ʷʰˤ[e]rʔ until Chinese-internal evidence appears that the vowel was *i.

OC *-r : Bodo -y : Lushai -i : Karen *-i/-j : Tib. zero Table

8

shows that Bodo, Lushai and Karen have *-i or *-j corresponding to OC *-r, Tib. zero. We know that Tib. has lost final *-j. Thus in Tibetan, loss of ST *-r in all likelihood occurred through a *-j stage.

"Words with OC *-r that were borrowed early into other languages sometimes lack the final -n that we find in Middle Chinese, suggesting either that the donor variety of Chinese was one that had changed *-r to *-j, or else that the borrowing language had no[r]. For example, for the cyclical sign 辰 chén < *[d]ər 'fifth earthly branch', Li (1945:336) cites these forms from Tai languages: Ahom shi, Lü si 1 , Dioi chi

. Khmu (a minority Mon-Khmer language spoken primarily in Laos) has /síi/(Damrong and Lindell 1994:104). If the original coda had been [n], there is no reason it could not have been borrowed with[n]."Baxter and Sagart (2014:255).

This form occurs as part of a disyllable ɕin 24 hɔ

'difficult, hard' < 辛苦 in a text at the end of Li's book(Li Jinfang 1999:259). Both syllables are consistent with the Early II layer. This word was probably borrowed as a disyllabic morpheme.

OC *-eŋ is reflected as /i/: tʰi 55 'hear' < *l^ˤeŋ (tone !); ɕi³³ 'star' < *s-tsʰˤeŋ; mi31 'name'< *C.meŋ.

The first example is doubtful (although the comparison itself is valid) because MC mwon is equally compatible with OC *mˤər, which is a good match for Tib. dmar. We should amend our Old Chinese reconstruction from *mˤur to *mˤ[ə]r. The second example is in all likelihood erroneous (the comparison itself is erroneous) because 君 jūn 'lord; ruler' relates to the action of governing (尹 *m-qurʔ > ywinX > yǐn 'govern; governor') rather than to honor. Thus there are no secure examples of OC *-ir, *-ur : corresponding to Tib. words ending in *-r.

Complementary distribution

The two correspondences, Hill's OC *-r : Tib. -r and my OC *-r : Tib. zero are thus in


Anna Marciniak-Czochra

Thomas Stiehl

Renormalization group second order approximation for singularly perturbed nonlinear ordinary differential equations

We consider a two time scale nonlinear system of ordinary differential equations. The small parameter of the system is the ratio ε of the scales. We search for an approximation involving only the slow time unknowns and valid uniformly for all times at order O(ε 2 ). It is a classical problem, studied using the Tikhonov's singular perturbation theorem. We develop an approach leading to a higher order approximation using the renormalization group (RG) method. We apply it in two steps. In the first step we show that the RG method allows to approximate the fast time variables by their RG expansion taken at the slow time unknowns. Next we study the slow time equations, where the fast time unknowns are replaced by their RG expansion and show the second order uniform error estimate. The procedure is computationally less demanding than the classical Vasil'eva-O'Malley expansion and allows a higher order extension of Hoppensteadt's result on the Tikhonov singular perturbation theorem for infinite times.

Introduction

In 1994, Chen, Goldenfeld and Oono proposed a simple and unified method for treating asymptotic perturbation problems [START_REF] Chen | Renormalization group theory for global asymptotic analysis[END_REF], [START_REF] Chen | Renormalization group and singular perturbations: Multiple scales, boundary layers, and reductive perturbation theory[END_REF]. This procedure originates from the quantum field theory and has been called the renormalization group (RG) method. In references [START_REF] Chen | Renormalization group theory for global asymptotic analysis[END_REF] and [START_REF] Chen | Renormalization group and singular perturbations: Multiple scales, boundary layers, and reductive perturbation theory[END_REF] it is demonstrated that the RG method has several advantages compared with the conventional methods of multiscale expansion. One advantage of the RG method is that its starting point is a naive perturbation expansion. Therefore, special intuition or experience to guess the asymptotic expansion are not needed. The RG method applies to solving linear and nonlinear boundary layer problems. It is well-known that straightforward perturbation expansions lead to the appearance of secular terms, i.e., terms exhibiting unbounded growth in time. Using the RG method, these terms are renormalized and lead to the RG equation which yields correct asymptotic expansions. Unexpectedly, the RG method seems to automatically identify all scales relevant to the problem [START_REF] Chen | Renormalization group theory for global asymptotic analysis[END_REF], [START_REF] Chen | Renormalization group and singular perturbations: Multiple scales, boundary layers, and reductive perturbation theory[END_REF]. Another remarkable property of the RG method is that it is sufficient to construct only the inner expansion. The RG method introduced in Ref. [START_REF] Chen | Renormalization group theory for global asymptotic analysis[END_REF], [START_REF] Chen | Renormalization group and singular perturbations: Multiple scales, boundary layers, and reductive perturbation theory[END_REF] has been further developed in Ref. [18], [START_REF] Mudavanhu | A New Renormalization Method for the Asymptotic Solution of Weakly Nonlinear Vector Systems[END_REF] and studied rigorously in Ref. [START_REF] Lee Deville | Analysis of a renormalization group method and normal form theory for perturbed ordinary differential equations[END_REF] and in a series of articles [START_REF] Chiba | C 1 approximation of vector fields based on the renormalization group method[END_REF]- [START_REF] Chiba | Extension and Unification of Singular Perturbation Methods for ODEs Based on the Renormalization Group Method[END_REF] by Chiba. The mathematically rigorous results apply to systems of ordinary differential equations of the form

dx dt = F x + εh(x, t),
where ε > 0 is a small parameter and F is a matrix with purely imaginary eigenvalues. Treating this problem using the naive expansion approach results in secular terms and the obtained asymptotic expansions are not valid for time intervals of the length

T = O( 1 ε
). In contrast to the naive expansion approach, the RG method provides a good approximation also for long times as shown in [START_REF] Chiba | C 1 approximation of vector fields based on the renormalization group method[END_REF]- [START_REF] Lee Deville | Analysis of a renormalization group method and normal form theory for perturbed ordinary differential equations[END_REF].

We note a number of papers on the application of the RG method to partial differential equations and, in particular, to the geostrophic flows, see [START_REF] Moise | Renormalization Group Method. Applications to Partial Differential Equations[END_REF], [START_REF] Temam | On the solutions of the renormalized equations at all orders[END_REF], [START_REF] Ziane | On a certain renormalization group method[END_REF]. For a study of the shadow limit for systems of ODEs coupled with nonlinear parabolic equations, we refer to [START_REF] Marciniak-Czochra | Shadow limits via the renormalization group method and the center manifold method[END_REF]. In the remainder of this paper we focus on systems of nonlinear ODEs.

We consider the asymptotic analysis of the ODE system with two characteristic times. Their ratio is equal to a small parameter ε > 0. Motivated by biological applications, we focus on the Cauchy problem

du ε dt = f (u ε , v ε ), t > 0, u ε (0) = u 0 ; (1) Tikh1A ε dv ε dt = -αv ε + Φ(u ε , v ε ), t > 0, v ε (0) = v 0 .
(2) Tikh1B

The asymptotic analysis of problem (1)-( 2) attracted a considerable interest in the literature.

The convergence to the solutions of the reduced system -αv + Φ(u, v) = 0 and

du dt = f (u, v), (3) Tikh1BS
was established for the stable roots v = φ(u) in the seminal theorem by Tikhonov (see e.g., [START_REF] Banasiak | Methods of small parameter in mathematical biology[END_REF] and [START_REF] Tikhonov | Translation of: Differentsial'nye uravneniiya)[END_REF]). It holds for finite time intervals and its extension to the infinite time interval is due to Hoppensteadt [START_REF] Hoppensteadt | Singular perturbations on the infinite interval[END_REF]. In applications, it is important to derive an error estimate for the approximation. For finite time intervals error estimates at all orders are achieved through boundary layer corrections by Vasil'eva and colleges (see [START_REF] Vasil'eva | The Boundary Function Method for Singular Perturbed Problems[END_REF]). The work of Hoppensteadt [START_REF] Hoppensteadt | Properties of Solutions of Ordinary Differential Equations with Small Parameters[END_REF] presents analogous estimates and also formulates, without proof, results for infinite time intervals. For a complete discussion of the subject we refer to a book by Hoppensteadt [START_REF] Hoppensteadt | Analysis and Simulation of Chaotic Systems[END_REF].

As observed by Chen, Goldenfeld and Oono in [START_REF] Chen | Renormalization group and singular perturbations: Multiple scales, boundary layers, and reductive perturbation theory[END_REF], the approach of Hoppensteadt [START_REF] Hoppensteadt | Properties of Solutions of Ordinary Differential Equations with Small Parameters[END_REF], O'Malley [START_REF] O'malley | Singular perturbation methods for ordinary differential equations[END_REF] and Vasil'eva [START_REF] Vasil'eva | The Boundary Function Method for Singular Perturbed Problems[END_REF] requires simultaneous expansion in all equations and leads to complex and non-straightforward calculations. Unlike the classical approach, the RG method allows to apply the asymptotic expansion first to equation (2) and then, independently, to equation [START_REF] Agarwal | An Introduction to Ordinary Differential Equations[END_REF]. After renormalization of the initial condition, a Cauchy problem for a single system of ODEs is obtained. With such choice of approximation, the RG method provides a uniform O(ε 2 ) error estimate on [T, +∞], T = O [START_REF] Agarwal | An Introduction to Ordinary Differential Equations[END_REF]. In addition, we construct effective initial conditions leading to a global O(ε 2 ) estimate.

Another powerful technique for studying system (1)-( 2) is the geometric singular perturbation theory (GSPT). The GSPT requires using normal forms. It applies to compact normally hyperbolic submanifolds of F (u, v) = -αv + Φ(u, v) = 0 , i.e., submanifolds consisting of roots of the nonlinear function F in which the eigenvalues of the gradient of F have non-zero real parts. For a detailed presentation of the method we refer to [START_REF] Kaper | An introduction to geometric methods and dynamical systems theory for singular perturbation problems. analyzing multiscale phenomena using singular perturbation methods[END_REF] and the references therein. The method has similarities to the Tikhonov theory but it is more general. The connection between RG and GSPT is discussed in [START_REF] Lee Deville | Analysis of a renormalization group method and normal form theory for perturbed ordinary differential equations[END_REF]. It is shown that RG provides a systematic method for finding normal forms for large classes of finite-dimensional vector fields. It has advantages over the normal form theory, since it allows systematic identification of secular terms by inspection of naive asymptotic expansions.

In this paper, using the RG theory, we derive an O(ε 2 ) approximation of ( 1)-( 2) that involves only the variables of the slow equations [START_REF] Agarwal | An Introduction to Ordinary Differential Equations[END_REF]. We obtain a Tikhonov type quasi-stationary approximation of order ε 2 for all times. Compared with the seminal work by Hoppensteadt [START_REF] Hoppensteadt | Properties of Solutions of Ordinary Differential Equations with Small Parameters[END_REF] and his monograph [START_REF] Hoppensteadt | Analysis and Simulation of Chaotic Systems[END_REF] our results are new, since we provide a proof of the error estimate of order ε 2 on the infinite interval. Another advantage of our construction is in its simplicity compared to the boundary layers calculations in the classical approach of Vasil'eva. Our presentation is self-contained and accessible to non-specialists in dynamical systems.

The paper is organized as follows. We formally derive the RG equation in Section 2. In Section 3 we derive the effective equations of problem ( 1)-( 2) using the center manifold approach. We notice that they are identical to the RG approximation from Section 2. Next, in Section 4 we introduce the setting in which the reduction to ( 3) is achieved at the order O(ε 2 ) for all times. Our construction works for any isolated stable root v of -αv + Φ(u, v) = 0. In Section 5 we present examples from the biosciences. The technical part of our results is postponed until Section 6.

RG approach to the singular perturbation

Sec2

We study Cauchy's problem (1)- [START_REF] Banasiak | Methods of small parameter in mathematical biology[END_REF]. The nonlinearities f and Φ are defined on R m+1 , m ≥ 1, and take values respectively in R m and R 1 . It is assumed that they are C 2 with bounded derivatives and that problem (1)-(2) has a unique globally defined smooth solution. α is a positive constant.

In order to apply the renormalization group (RG) approach, we change the time scale by setting τ = t/ε. System (1)-( 2) becomes

du ε dτ = εf (u ε , v ε ), dv ε dτ = -αv ε + Φ(u ε , v ε ), u ε (0) = u 0 , v ε (0) = v 0 . ( 4
) TikhRG1
We now proceed by the usual RG method steps.

1. We pretend that the problem can be solved as a regular perturbation problem and try a naive expansion

u ε (τ ) = u 0 (τ ) + εu 1 (τ ) + ε 2 u 2 (τ ) + . . . , v ε (τ ) = v 0 (τ ) + εv 1 (τ ) + . . . It gives d dτ (u 0 (τ ) + εu 1 (τ ) + ε 2 u 2 (τ ) + . . . ) = εf (u 0 , v 0 )+ ε 2 (∇ u f (u 0 , v 0 )u 1 (τ ) + ∂ v f (u 0 , v 0 )v 1 (τ )) + . . . , d dτ (v 0 (τ ) + εv 1 (τ ) + . . . ) = -α(v 0 (τ ) + εv 1 (τ ) + . . . ) + Φ(u 0 , v 0 )+ ε{∇ u Φ(u 0 , v 0 )u 1 (τ ) + ∂Φ ∂v (u 0 , v 0 )v 1 (τ )} + . . . .
We sort terms with respect to the powers of ε. Therefore, at order zero we have

d dτ v 0 (τ ) = -αv 0 + Φ(u 0 (τ ), v 0 (τ )) and d dτ u 0 (τ ) = 0. ( 5
) RGSP1
We set u 0 (0) = A and v 0 (0

) = B (6) RGSP2 Then u 0 (τ ) = A. ( 7
) RGSP2A
In order to continue we have to make additional assumptions on Φ.

First, we restrict ourselves to the case with -αv + Φ(A, v) = 0 having a unique solution. This assumption follows the Tikhonov theory and we make it here to streamline the presentation. The assumption will be relaxed in Section 4.

Assumption 1 It holds

-κ 0 ≤ -ξ := -α + ∂Φ ∂v ≤ -κ < 0 (8) HRG
in R m+1 for some positive constants κ 0 and κ 1 .

We denote by φ(A) the unique root of -αv + Φ(A, v) = 0.

Dec1

Lemma 1 Under hypothesis [START_REF] Chiba | Analysis of a renormalization group method and normal form theory for perturbed ordinary differential equations[END_REF] we have the estimate Proof. Denote by sign µ (•), µ ∈ R regularizations of the sign function that converge to the sign function for µ passing to zero. Multiplying equation ( 5) by regularized sign (v 0 (τ ) -φ(A)) and passing to the zero limit of the regularization parameter µ, we obtain the estimate [START_REF] Chiba | Extension and Unification of Singular Perturbation Methods for ODEs Based on the Renormalization Group Method[END_REF].

|v 0 (τ ) -φ(A)| ≤ Ce -κτ , ( 9
We now consider terms of order O(ε) and obtain

d dτ v 1 (τ ) = -αv 1 (τ ) + ∇ u Φ(A, v 0 )u 1 (τ ) + ∂Φ ∂v (A, v 0 )v 1 (τ ), ( 10
) RGSP3 d dτ u 1 (τ ) = f (A, v 0 (τ )). ( 11
) RGSP4
In analogy to Lemma 1, we have:

Dec11

Lemma 2 Under hypothesis [START_REF] Chiba | Analysis of a renormalization group method and normal form theory for perturbed ordinary differential equations[END_REF] there is a constant C 1u such that the following estimate holds

|u 1 (τ ) -τ f (A, φ(A)) -C 1u | + | d dτ (u 1 (τ ) -τ f (A, φ(A)))| ≤ Ce -κτ . (12) Decayu1 Consequently, u 1 (τ ) = τ f (A, φ(A)) + C 1u + O(e -κτ ).
Proof. In analogy to Lemma 1.

It remains to calculate v 1 and u 2 . We use equation [START_REF] Lee Deville | Analysis of a renormalization group method and normal form theory for perturbed ordinary differential equations[END_REF] to calculate v 1 . Using ( 9) and ( 12), we find out that equation [START_REF] Lee Deville | Analysis of a renormalization group method and normal form theory for perturbed ordinary differential equations[END_REF] can be written as

d dτ v 1 (τ ) = -ξ A v 1 (τ ) + ∇ u Φ(A, φ(A))(τ f (A, φ(A)) + C 1u ) + F 1 , ( 13
)
RGSP3AA
where

ξ A = α - ∂Φ ∂v (A, φ(A)) and |F 1 (τ )| ≤ Ce -κ 1 τ , κ 1 < κ.
Note that ξ A is strictly positive. Therefore, we have

v 1 (τ ) = v 1 (0) exp{-ξ A τ } + ∇ u Φ(A, φ(A)) • ∫ τ 0 exp{-ξ A (τ -y)}(C 1u +yf (A, φ(A))) dy + ∫ τ 0 exp{-ξ A (τ -y)}F 1 (y) dy = ∇ u Φ(A, φ(A)) • 1 ξ A ( τ f (A, φ(A)) + C 1u - f (A, φ(A)) ξ A ) + O(e -κ 1 τ ), κ 1 < κ. ( 14
) U10
At order O(ε 2 ) for u 2 we have

d dτ u 2 (τ ) = ∇ u f (u 0 , v 0 )u 1 (τ ) + ∂ v f (u 0 , v 0 )v 1 (τ ) (15) RGSP4A
Using (9), ( 12) and ( 14), we find out that equation ( 15) can be written as

d dτ u 2 (τ ) = ∇ u f (A, φ(A))(τ f (A, φ(A)) + C 1u ) + ∂ v f (A, φ(A)) ( +∇ u Φ(A, φ(A)) • 1 ξ A ( τ f (A, φ(A))+C 1u - f (A, φ(A)) ξ A ) ) + F 2 , ( 16
) RGu2
where

|F 2 (τ )| ≤ Ce -κ 1 τ , κ 1 < κ. Therefore, we have u 2 (τ ) = τ 2 2 ( ∇ u f (A, φ(A))f (A, φ(A))+ ∂ v f (A, φ(A))∇ u Φ(A, φ(A)) • 1 ξ A f (A, φ(A)) ) + τ ( ∇ u f (A, φ(A))C 1u + ∂ v f (A, φ(A))∇ u Φ(A, φ(A)) • 1 ξ A ( C 1u - f (A, φ(A)) ξ A ) ) + C 2u + O(e -κ 1 τ ), κ 1 < κ, (17) U10u2
where C 2u is an arbitrary constant vector.

The approximation is now

u ε (τ ) = u 0 (τ ) + εu 1 (τ ) + ε 2 u 2 (τ ) + O(ε 3 ), ( 18
) RGSP5A v ε (τ ) = v 0 (τ ) + εv 1 (τ ) + O(ε 2 ). ( 19
) RGSP5B
In most case, we detect immediately the presence of secular terms in u 1 . They can also appear in u 2 and v 1 .

2. The idea of the normalization is to introduce an arbitrary time µ, split τ as τ -µ + µ and absorb the terms containing µ into the renormalized counterpart A(µ) of A. We introduce the renormalization constant

Z 1 = 1 + a 1 ε + a 2 ε 2 . The coefficient A is replaced by (1+a 1k ε+a 2k ε 2 )A k (µ), k = 1, . . . , m.
The coefficients a 2 and a 1 are chosen to eliminate the terms containing µ at order O(ε 2 ).

Our approximation ( 18)-( 19) becomes

u k (τ ) = u 0k (τ ) + εu 1k (τ ) + ε 2 u 2k (τ ) = (1 + a 1k ε + a 2k ε 2 )A k (µ)+ ε(τ -µ + µ)f k ( A(µ)(1 + a 1 ε)), φ(A(µ)(1 + a 1 ε)) ) + εC 1u,k + ε 2 u 2k (τ ), ( 20
) RGSP6A v = v 0 (µ) + εv 1 (µ) = φ(A(µ)(1 + a 1 ε))+ ε∇ u Φ(A, φ(A)) • 1 ξ A ( (τ -µ + µ)f (A, φ(A)) + C 1u - f (A, φ(A)) ξ A ) . ( 21
) RGSP6B
We choose a 1k such that the term µf k (Z 1 A, φ(Z 1 A)) is eliminated at the leading order. Hence we have

a 1k (µ)A k (µ) + µf k (A, φ(A)) + C 1u,k = 0, implying a 1k = - µ A k f k (A, φ(A)) - C 1u,k A k , k = 1, . . . , m. ( 22
) RGSP8
Therefore, for u 1 we have the following expression

u 1 (τ ) = (τ -µ)f ((1 + a 1 ε))A, φ((1 + a 1 ε))A)) + R = (τ -µ)f (A, φ(A)) -ε(τ -µ)(∇ u f (A, φ(A))C 1u + ∂ v f (A, φ(A))∇ A φ(A) • C 1u )- ε(τ -µ)µ ( ∇ u f (A, φ(A))f (A, φ(A)) + ∂ v f (A, φ(A)) ξ A ∇ u Φ(A, φ(A)) • f (A, φ(A)) ) +R, (23) RGSP8A
where R and its derivative with respect to τ are exponentially small in τ .

Next, expressions ( 14) and ( 22) yield

v 0 (µ) + εv 1 (µ) = φ(A) + ε ( ∇ A φ(A) • a 1 A + + ∇ u Φ(A, φ(A)) • 1 ξ A ( (τ -µ + µ)f (A, φ(A)) + C 1u - f (A, φ(A)) ξ A ) ) = φ(A) + ε(τ -µ - 1 ξ A )∇ u Φ(A, φ(A)) • f (A, φ(A)) ξ A . ( 24
) U1
We now transform τ into τ -µ+µ and τ 2 into τ 2 -µ 2 +µ 2 in the expression (17) for u 2 (µ).

Inserting formulas [START_REF] Nakata | Stability analysis of multicompartment models for cell production systems[END_REF] and renormalized ( 17) into (20) yields u k (τ ) without secular terms in µ. The terms of order O(ε 2 ), containing only µ and µ 2 , are to be eliminated and only terms containing τ -µ and τ 2 -µ 2 remain. We achieve this goal by choosing appropriate a 2k . After recalling that

∇φ(A)(α - ∂Φ ∂v (A, φ(A)) = ∇ u Φ(A, φ(A)),
we obtain the following expression for u(τ ):

u(τ ) = A(µ) + ε(τ -µ)f (A, φ(A)) -ε 2 (τ -µ) ∇ u Φ(A, φ(A)) • f (A, φ(A)) (ξ A ) 2 ∂ v f (A, φ(A)) + ε 2 2 (τ -µ) 2 ( ∇ u f (A, φ(A))f (A, φ(A))+ ∂ v f (A, φ(A)) ξ A ∇ u Φ(A, φ(A)) • f (A, φ(A)) ) . ( 25
) fullrenor 3.
The parameter µ is artificial and the solution does not depend on it. Therefore, it has to hold ∂u(τ ) ∂µ | µ=τ = 0 for all τ . After noticing that terms multiplying ε 2 µ cancel, ∂u(τ ) ∂µ | µ=τ = 0 implies the RG equation

dA dτ = εf (A, φ(A)) -ε 2 ∇ u Φ(A, φ(A)) • f (A, φ(A)) (α - ∂Φ ∂v (A, φ(A))) 2 ∂ v f (A, φ(A)). ( 26
) RGeqSP1
Returning to the original variable t = ετ , we obtain

dA dt = f (A, φ(A)) -ε ∇ u Φ(A, φ(A)) • f (A, φ(A)) (α - ∂Φ ∂v (A, φ(A))) 2 ∂ v f (A, φ(A)), ( 27
) RGeqSP1A
which is the RG equation we were looking for. For t = O(1) the initial time layer effects became negligible and the approximation is expressed by A(t).

It remains to check the value of the derivative of v(τ

) = v 0 (µ) + εv 1 (µ) with respect to µ at µ = τ . A direct calculation yields ∂v(τ ) ∂µ | µ=τ = ∇φ(A) dA dτ -ε f (A, φ(A)) • ∇ u Φ(A, φ(A)) α - ∂Φ ∂v (A, φ(A)) - ε d dµ ( f (A, φ(A)) • ∇ u Φ(A, φ(A)) (α - ∂Φ ∂v (A, φ(A))) 2 ) = O(ε 2 ).
Therefore, the choice

u ε (t) ≈ A RG (t), ( 28
) RGAppr1 v ε (t) ≈ φ(A RG (t)) -ε f (A RG , φ(A RG )) • ∇ u Φ(A RG , φ(A RG )) (αA RG -Φ(A RG , φ(A RG ))) 2 , ( 29
) RGAppr2
where A RG is expressed by [START_REF] Vasil'eva | The Boundary Function Method for Singular Perturbed Problems[END_REF], provides the requested O(ε 2 ) approximation, for times O(1) ≤ t ≤ O(1/ε). The correct behavior for small times is described by the initial time layers, as in the classical literature (see e.g., [START_REF] O'malley | Singular perturbation methods for ordinary differential equations[END_REF], [START_REF] Tikhonov | Translation of: Differentsial'nye uravneniiya)[END_REF] and [START_REF] Vasil'eva | The Boundary Function Method for Singular Perturbed Problems[END_REF]).

Center Manifold Theorem approach to obtain the effective model

SecCMF

In this section we present a formal derivation of the effective equations using the Center Manifold Theorem.

The usual way to study the qualitative behavior of the solutions of the equation

dy dt = Ψ(y), Ψ(0) = 0, ( 30
) CMT1
is to calculate the Jacobian matrix A = ∇Ψ| y=0 and to study the linearized problem

dy dt = Ay + G(y), G(0) = 0, ∇G(0) = 0. ( 31
) CMT2
If all the eigenvalues of the matrix A = ∇Ψ| y=0 have non-zero real part, then the Hartman-Grobman theorem applies and the behavior of equation ( 30) in a small neighborhood of 0 is described by equation (31).

If there exist eigenvalues with zero real part, a center manifold W c is spanned by the corresponding (generalized) eigenvectors. In this case the Center Manifold Theorem is used to restrict the analysis to such a "center manifold". For an introduction to the center manifold technique, we quote the classical book by Carr [START_REF] Carr | Applications of Centre Manifold Theory[END_REF].

Our goal is to study behavior of the solution to problem (1)-( 2), namely

du ε dt = f (u ε , v ε ), t > 0, u ε (0) = u 0 ; dv ε dt = - α ε v ε + 1 ε Φ(u ε , v ε ), t > 0, v ε (0) = v 0 ,
for ε → 0. We notice immediately the large difference between the modulus of the eigenvalues related to equation [START_REF] Banasiak | Methods of small parameter in mathematical biology[END_REF] and those related to equation [START_REF] Agarwal | An Introduction to Ordinary Differential Equations[END_REF]. Therefore, it is to be expected that equation (1) will play the role of the center manifold for problem (1)- [START_REF] Banasiak | Methods of small parameter in mathematical biology[END_REF].

More precisely, the solution for (2) is given by

v ε - 1 α Φ(u ε , v ε ) = v 0 e -αt/ε + ( 1 ε e -αt/ε - 1 α δ(t)) * Φ(u ε , v ε ), (32) Eqslave1
where * denotes the convolution product in time and δ(t) the delta distribution centered at t. Under mild assumptions on Φ, the right hand side converges to zero, as ε → 0. We conclude

that v ε ≈ 1 α Φ(u ε , v ε )
, for small ε. Hence, it is reasonable to assume that v lies on an invariant manifold, which we represent as v = Z(u ε ). On the invariant manifold it holds

dv dt = d dt Z(u) = ∇Z d dt u, ∀t.
We substitute the above identity into (2) and obtain

-αZ(u) + Φ(u, Z(u)) = ε∇Z • f (u, Z(u)). ( 33
) CMeq1
We expand Z as

Z = Z 0 + εZ 1 + O(ε 2 ).
Sorting terms with respect to their orders of ε, we obtain

O(1) : -αZ 0 + Φ(u, Z 0 ) = 0; O(ε) : -αZ 1 + ∂ v Φ(u, Z 0 )Z 1 = ∇Z 0 • f (u, Z 0 ).
The first equation implies that Z 0 is a root of v -1 α Φ(u, v) = 0. The second equation implies that

Z 1 = ∇Z 0 • f (u, Z 0 ) -α + ∂ v Φ(u, Z 0 ) . Hence v = Z 0 + ε ∇Z 0 • f (u, Z 0 ) -α + ∂ v Φ(u, Z 0 ) . ( 34
) Approv
To obtain the ODE for u, we substitute expression (34) into (1):

du dt = f ( u, φ(u) -ε ∇ u Φ(u, φ(u)) • f (u, φ(u)) (α - ∂Φ ∂v (u, φ(u))) 2 ) = f (u, φ(u)) -ε ∇ u Φ(u, φ(u)) • f (u, φ(u)) (α - ∂Φ ∂v (u, φ(u))) 2 ∂ v f (u, φ(u)) + O(ε 2 ). ( 35
) RGeqSP1AC
We note that equation ( 35) is identical to equation [START_REF] Vasil'eva | The Boundary Function Method for Singular Perturbed Problems[END_REF].

Remark 1 It is possible to transform the above intuitive argument into a rigorous one. Nevertheless, most results obtained using the Center Manifold Theorem describe dynamics of small solutions around the origin. In our problem we would like to prove that for small values of the parameter ε we remain for all times in an ε k -neighborhood of the solution of equation [START_REF] Vasil'eva | The Boundary Function Method for Singular Perturbed Problems[END_REF].

The reduction to a neighborhood of zero is not straightforward and we prefer presenting a direct proof, based on the computations from Section 2 and the boundary layers, adopted from Tikhonov's theory.

Main results

Sec4

We start by making basic assumptions on our singularly perturbed system.

Ass1 Assumption 2 Let (f , Φ) ∈ C 3 (R m+1 ) m+1 , m ≥ 1,

and let the algebraic equation

-αy + Φ(x, y) = 0 in R, (36) Root1
have at least one isolated real root y = y(x).

Ass2 Assumption 3 We suppose that the reduced problem

-αv + Φ(A, v) = 0; dA dt = f (A, v), t > 0, and
A(0) = u 0 , ( 37
) Reduced1
has a smooth bounded solution {A, v} on R + and that v is an isolated real root.

Ass3 Assumption 4 We suppose that the chosen root v satisfies

ξ = α -∂ y Φ(x, y)| (x=A,y=v) ≥ M 0 > 0 on R + , (38) Root2
i.e., that v is a stable root.

In addition, we suppose that the initial datum v 0 is in the basin of attraction of the root v from (37) (i.e., for all δ > 0, there exists t(δ) > 0 such that the solution w of the initial value problem d dt w = -αw + Φ(A(t), w), w(0) = v 0 satisfies |w(t) -φ(A(t))| ≤ δ for t > t(δ)). We furthermore assume that on the interval

I = [min{φ(u 0 ), v 0 }, max{φ(u 0 ), v 0 }] we have α -∂ y Φ(x, y)| (x=u 0 ,y∈I) ≥ M 0 2 > 0. ( 39
) Root2A
Next, we introduce the initial layer ζ 0 by

dζ 0 dτ = -αζ 0 + Φ(u 0 , φ(u 0 ) + ζ 0 (τ )) -Φ(u 0 , φ(u 0 )), ζ 0 (0) = v 0 -φ(u 0 ) (40) Inn2pbz
(to be compared with ( 9)). Note that by Assumption 4, φ(u 0 ) is uniquely defined and applying Lemma 1 implies

|ζ 0 (τ )| ≤ Ce -M 0 τ . ( 41
) Dec2
Note that ζ 0 (t/ε) is the classical initial layer function from [START_REF] Hoppensteadt | Properties of Solutions of Ordinary Differential Equations with Small Parameters[END_REF], [START_REF] O'malley | Singular perturbation methods for ordinary differential equations[END_REF] and [START_REF] Vasil'eva | The Boundary Function Method for Singular Perturbed Problems[END_REF].

For an error estimate another assumption is needed:

Ass4 Assumption 5 We suppose that the matrix

( ∇ u f (x, φ(x)) + ∂ v f (x, φ(x)) ⊗ ∇φ(x) ) | x=A (42) Stabil1
satisfies the exponential dichotomy assumption on R + , i.e., we assume that a fundamental solution S of the linear system corresponding to (42) fulfills S(0) = I, ||S(t)S -1 (τ )|| ≤ Ke -κ(t-s) for s ≤ t ≤ ∞ and appropriate positive constants K, κ.

Dichotomy Remark 2 We recall that this is a special case of the exponential dichotomy assumption, stating that there exists a projection P , such that ||S(t)P S -1 (τ )|| ≤ Ke -κ(t-s) for s ≤ t ≤ ∞ and ||S(t)(I -P )S -1 (τ )|| ≤ Le µ(t-s) for s ≥ t ≥ -∞ for appropriate positive constants K, L, κ, µ, [START_REF] Hoppensteadt | Properties of Solutions of Ordinary Differential Equations with Small Parameters[END_REF], [START_REF] Lin | Linear Systems Exponential Dichotomy and Structure of Sets of Hyperbolic Points[END_REF].

Remark2 Remark 3 Let b > 0 and

T b,λ = Int ∪ t∈R + {A(t) + x, φ(A) + y} {|x j |<b, j=1,...,m, y∈I λ (t)} , (43) tube
where λ(t) is a smooth curve, defined on R + and such that λ(0) = |φ(u 0 )-v 0 | and λ(t) = λ 0 > 0 for t ≥ t 0 . The set I λ is given by

I λ (t) = { (-b, λ(t)), for v 0 ≥ φ(u 0 ), (-λ(t), b), for v 0 < φ(u 0 ). ( 44
) ilam
Since A is bounded and continuous, and because of (39) in Assumption 4, there exists a tubular domain T b,λ defined by (43) and depicted in Fig. 1, such that (38) and (42) hold true in T b,λ . Note that the existence of exponential dichotomies of the considered form is preserved under small continuous perturbations or continuous perturbations vanishing as t → ∞, see [START_REF] Lin | Linear Systems Exponential Dichotomy and Structure of Sets of Hyperbolic Points[END_REF] or [START_REF] Agarwal | An Introduction to Ordinary Differential Equations[END_REF] respectively. Note that for all times |ζ 0 (t/ε)| ≤ λ(t). We now define C ∞ -cutoff extensions of the nonlinearities f and Φ. Let ξ ∈ R and g a C ∞ -function such that

g(ξ) =    ξ,
for |ξ| ≤ 1/2; 1, for ξ ≥ 1; -1, for ξ ≤ -1.

(45) smoothy

Then we replace Φ by Φ, where

Φ(x 1 , . . . , x m+1 , t) = Φ(x 1 , . . . , xm+1 ), (x 1 , . . . , x m+1 , t) ∈ R m+1 × R + , ( 46
) Extphi and xj = g( x j -A j (t) b )b + A j (t), (47) Extphi1
for j = 1, . .

. , m and xm+1 = g( x m+1 -φ(A(t)) λ(t) )λ(t) + φ(A(t)). (48) Extphi2

The function Φ is smooth, globally Lipschitz and satisfies condition (38) for (x, y) ∈ R m+1 and t ∈ R + .

We replace f by f , where

f (x 1 , . . . , x m+1 , t) = f (x 1 , . . . , xm+1 ), (x 1 , . . . , x m+1 , t) ∈ R m+1 × R + , ( 49
du ε dt = f (u ε , v ε , t), t > 0, u ε (0) = u 0 ; (50) Tikh1E ε dv ε dt = -αv ε + Φ(u ε , v ε , t), t > 0, v ε (0) = v 0 . (51) Tikh1F
has a unique C 1,1 solution on [0, T ], for all 0 < T < +∞. Furthermore, the C 1,1 norm is uniformly bounded.

Remark 4

The above C ∞ -cutoff extension of the nonlinearities is motivated by the similar extensions in the center manifold theory. For more details we refer to the seminal book of Carr [START_REF] Carr | Applications of Centre Manifold Theory[END_REF], p.16-17.

We study the error function for equation (51) and define

V ε (t): = v ε (t) -φ(u ε (t), t) -ζ0 ( t ε ), t > 0, (52) V1
where φ is the modification of the root φ of equation (36), obtained after replacing Φ by Φ and ζ0 is the solution of problem (40) obtained after replacing Φ by Φ.

Prop0 Proposition 1 |V ε (t)| ≤ Cε, 0 ≤ t, ( 53
) EstV0
where C does not depend on the extension.

Proof. See Section 6.2.

Next, we study the error function for u ε and define

δ ε (t): = u ε (t) -A(t), t ≥ 0. ( 54
) Erru1
Under the above assumptions, the error function δ ε satisfies the following Cauchy problem

dδ ε dt = f ( A + δ ε , φ(A + δ ε , t) + ζ 0 ( t ε ) + V ε , t ) -f (A, φ(A), t), t > 0; (55) Errorder1 δ ε (0) = 0. ( 56
) Tikh1G
The following result holds:

Prop2 Proposition 2 Problem (55)-( 56) has a unique solution δ ε satisfying

|δ ε (t)| ≤ Cε, 0 ≤ t, (57) Erru2
where C does not depend on the extension and the time interval, furthermore C ≤ λ(0).

Proof. See Section 6.2.

Cor1 Corollary 1 For ε ≤ ε 0 it holds f = f and Φ = Φ in a neighborhood of the solution. Hence the solutions of problem (50)-(51) coincide with the solution of problem ( 1)-( 2). They exist for all times.

Remark 5 In [START_REF] Hoppensteadt | Singular perturbations on the infinite interval[END_REF] the proof that the solution of problem ( 1)-( 2) remains in the tubular neighborhood T b,λ is carried out using a local Lyapunov functional. In our setting, this property results from the error estimate.

longtime Remark 6 Let us replace the exponential dichotomy from Assumption 5 by the assumption that the matrix

( ∇ u f (x, φ(x))+∂ v f (x, φ(x))⊗∇φ(x)
)

| x=A has a symmetric part with non-positive eigenvalues.

(58) Stabil1L Then for every γ > 0,

|δ ε (t)| ≤ Cε 1-γ , 0 ≤ t ≤ Cε -γ , ( 59
) Erru2T
where C does not depend on the extension and the time interval.

Proof. See Section 6.2.

Th1 Theorem 1 Under Assumptions 2-5, it holds

sup 0≤t≤+∞ |u ε (t) -A(t)| ≤ Cε, (60) ESTO1 |v ε (t) -φ(A(t))| ≤ C(ε + e -M 0 t/ε ), ∀t ∈ R + . ( 61
) ESTO2
Proof. The statement follows from Propositions 1 and 2.

Prop3 Proposition 3 Let

U ε = v ε -φ(u ε ) + ε α -∂Φ ∂v (u ε , φ(u ε )) d dt φ(u ε ) -ζ 0 (t/ε) -εζ 1 (t/ε),
where

ζ 1 fulfills d dt ζ 1 (t) = -αζ 1 (t), ζ 1 (0) = ∇φ(u 0 ) α -∂Φ ∂v (u 0 , φ(u 0 )) f (u 0 , v 0 ). ( 62
) zeta1
Under Assumptions 2-5, it holds

U ε (t) = U ε,1 (t) + U ε,2 (t), (63)
where U ε,1 are U ε,2 are differentiable functions fulfilling

||U ε,1 || L ∞ (R + ) ≤ Cε 2 , ||U ε,2 || L 1 (R + ) ≤ Cε 2 . ( 64
)
Proof. See Section 6.2.

Prop4 Proposition 4 Let Assumptions 2-5 be satisfied. Define δε (t) = u ε (t) -A(t) -εu 1 (t) -εψ 1 (t/ε) with

{ d dτ ψ 1 = f ( u 0 , φ(u 0 ) + ζ 0 (τ )) -f (u 0 , φ(u 0 )) ψ 1 (0) = - ∫ ∞ 0 ( f ( u 0 , φ(u 0 ) + ζ 0 (τ )) -f (u 0 , φ(u 0 )) ) dτ ( 65
)
and

         d dt u 1 = ( ∇ u f (x, y) + ∂ v f (x, y) ⊗ ∇φ(x) ) | (x=A,y=φ(A)) u 1 -∂ v f (A, φ(A)) ∇φ(A)•f (A,φ(A)) α-∂Φ ∂v (A,φ(A)) u 1 (0) = -ψ 1 (0). ( 66
)
Then, the following estimate holds for all t ≥ 0:

| δε (t)| ≤ Cε 2 .
Proof. See Section 6.2.

d dt A 1 = ( 2a 1 + kA 2 -1)pA 1 , A 1 (0) = u 0 1 (73) SC10 d dt A 2 = 2(1 - a 1 + kA 2 )pA 1 -dA 2 , A 2 (0) = u 0 2 (74) SC20 v = 1 1 + kA 2
has been extensively studied in references [START_REF] Getto | Global dynamics of twocompartment models for cell production systems with regulatory mechanisms[END_REF], [START_REF] Nakata | Stability analysis of multicompartment models for cell production systems[END_REF] and [START_REF] Stiehl | Characterization of stem cells using mathematical models of multistage cell lineages[END_REF], motivated by the different time scales inherent to the biological system. Our results provide a rigorous framework to show that in proximity of positive equilibria the distance between the solutions of the reduced system (73)-( 74) and the original system (70)-( 71) is bounded by Cε, and the distance of the solutions of the original system (70)-( 71) and the first order RG-approximation is bounded by Cε 2 for times t ≥ O(1).

ExampleLemma1 Lemma 3 Consider system (70)-( 72) and the corresponding 0th-order RG-approximation ( 73)-( 74). Both systems have a unique positive equilibrium. In the vicinity of this equilibrium Assumptions 2-5 are satisfied.

Proof. We have

α = 1, Φ(x, y) = -kx 2 y + 1, f (x, y) = [(2ay -1)px 1 , 2(1 -ay)px 1 -dx 2 ] T .
Assumption 2 is trivially fulfilled with the isolated root y(x) = φ(x) = 1 1+kx 2 . We note that A is nonnegative for nonnegative initial conditions. Assumption 4 is fulfilled, since α -

∂ y Φ(x, y)| x=A,y=φ(A) = 1 + kA 2 ≥ 1.
To show the boundedness of A, we consider q := A 1 A 2 , which fulfills the initial value problem q(0)

= u 0 1 /u 0 2 d dt q = (2 a 1+kA 2 -1)pq -2(1 -a 1+kA 2 )pq 2 + dq ≤ (2a -1)pq -2(1 -a)pq 2 + dq. Consequently, q > (2a-1)p+d 2 2(1-a)p =: Q implies d dt q < 0. There- fore, q ≤ max{Q, q(0)} =: K 1 . Then we obtain d dt A 1 ≤ ( 2a 1+kA 1 /K 1 -1)pA 1 . Therefore, A 1 > (2a-1)K 1 k implies d dt A 1 < 0. Consequently, A 1 ≤ max{u 0 1 , (2a-1)K 1 k } =: K 2 . This yields d dt A 2 ≤ 2pK 2 -dA 2 . Therefore, A 2 ≤ max{u 0 2 , 2pK 2 d }.
Local Lipschitz continuity and global boundedness imply existence of global solutions, Assumption 3 is thus satisfied. The positive equilibrium is given by ȳ =

1 1+kū ε,2 , Ā1 = ūε,1 = (2a-1)d kp , Ā2 = ūε,2 = 2a-1
k . The linearization of system (73)-(74) around the positive equilibrium has two eigenvalues with negative real parts, due to Vieta's theorem. Therefore, Assumption 5 is fulfilled for the constant solution A = Ā. Since Ā is locally stable it follows that solutions A of system (73)-( 74) with initial values in the vicinity of Ā converge exponentially to Ā. Due to stability of exponential dichotomies with respect to small continuous perturbations, Assumption 5 is fulfilled along trajectories of such solutions A.

The 1st-order RG approximation is given by the following Lemma.

Lemma 4 The 1st-order RG approximation of (70)-( 72) is given by

d dt A RG 1 = ( 2a 1 + kA RG 2 -1 ) pA RG 1 +2apA RG 1 k 2 ( 1 - a 1+kA RG 2 ) pA RG 1 -dA RG 2 ( 1 + kA RG 2 ) 3 ε (75) ARG1 d dt A RG 2 = [2 ( 1 - a 1 + kA RG 2 ) pA RG 1 -dA RG 2 ] ( 1 - 2apA RG 1 k (1 + kA RG 2 ) 3 ε ) (76) ARG2 A RG (0) = [u 0 1 , u 0 2 ] T + 2apu 0 1 ( v 0 -1 1+ku 0 2 ) 1 + ku 0 2 [1, -1] T ε. ( 77
) ARG3
Proof. We use the expressions for α, Φ, f , φ specified in the proof of Lemma 3, together with

φ(u) = 1 1+ku 2 , ∇φ(u) = [0, -k (1+ku 2 ) 2 ] T , ∂ v f (u, v) = 2apu 1 [1, -1] T . The ODEs for A RG 1 and A RG 2 follow from Theorem 2.
For the initial layer

ζ 0 we obtain ζ 0 (0) = v 0 -1 1+ku 0 2 , d dt ζ 0 = -(1 + ku 0 2 )ζ 0 and therefore ζ 0 (t) = ( v 0 -1 1+ku 0 2 ) e -(1+ku 0 2 )t . Furthermore, it holds f (u 0 , φ(u 0 ) + ζ 0 ) -f (u 0 , φ(u 0 )) = 2aζ 0 pu 0 1 [1, -1] T and ∫ ∞ 0 ( f (u 0 , φ(u 0 ) + ζ 0 ) -f (u 0 , φ(u 0 )) ) dt = ( v 0 -1 1+ku 0 2 ) 1 + ku 0 2 [2apu 0 1 , -2apu 0 1 ] T .
This yields the initial condition. Lemma 3 and Theorem 2 imply that the distance of solutions of the original system (70)-(71) and the first order RG-approximation (75)-( 77) is bounded by Cε 2 for times t ≥ O(1).

Proof of the main results

Sec5

Proof of Proposition 1

Proof. In order to simplify notation, we systematically skip tilde overbars in the proof.

By definition V ε (0) = 0 and

ε d dt V ε = -αv ε + Φ(u ε , v ε , t) -ε d dt φ(u ε , t) + αζ 0 (t/ε) + Φ(u 0 , φ(u 0 , t), t) -Φ(u 0 , φ(u 0 , t) + ζ 0 (t/ε), t) = -αV ε + Φ(u ε , v ε , t) -Φ(u ε , φ(u ε , t) + ζ 0 ( t ε ), t) -ε d dt φ(u ε , t) +α(φ(u 0 , t) + ζ 0 ( t ε )) + Φ(u ε , φ(u ε , t) + ζ 0 (t/ε), t) -α(φ(u ε , t) + ζ 0 ( t ε )) -Φ ( u 0 , φ(u 0 , t) + ζ 0 ( t ε ), t ) = -(α - ∂Φ ⋆ ∂v )V ε -ε d dt φ(u ε , t) + I(t) (78) mathcalVO with I(t) := Φ(u ε , φ(u ε , t) + ζ 0 (t/ε), t) -α(φ(u ε , t) + ζ 0 (t/ε)) -Φ ( u 0 , φ(u 0 , t) + ζ 0 (t/ε), t ) + α(φ(u 0 , t) + ζ 0 (t/ε)) (79) IT and ∂Φ ⋆ ∂v := ∂Φ(u ε , v ⋆ ) ∂v
, where v ⋆ is an intermediate value between v ε and φ(u ε ). In the above calculations we used equation αφ(x) = Φ(x, φ(x, t), t) for all x ∈ R m . Next, we use the smoothness of the extensions and estimate I(t) using the second order derivatives:

|I(t)| = ∫ 1 0 d dη [Φ(u ε (ηt), φ(u ε (ηt)) + ζ 0 (t/ε)) -α(φ(u ε (ηt)) + ζ 0 (t/ε))]dη = ∫ 1 0 t∇ u Φ(u, v)| u=uε(ηt),v=φ(uε(ηt))+ζ 0 (t/ε)) • d dτ u ε (τ )| τ =ηt -t∇ u Φ(u, v)| u=uε(ηt),v=φ(uε(ηt)) • d dτ u ε (τ )| τ =ηt +t∂ v Φ(u, v)| u=uε(ηt),v=φ(uε(ηt))+ζ 0 (t/ε)) ∇φ(u ε (ηt)) • d dτ u ε (τ )| τ =ηt -t∂ v Φ(u, v)| u=uε(ηt),v=φ(uε(ηt)) ∇φ(u ε (ηt)) • d dτ u ε (τ )| τ =ηt -tα∇φ(u ε (ηt)) • d dτ u ε (τ )| τ =ηt + tα∇φ(u ε (ηt)) • d dτ u ε (τ )| τ =ηt dη ≤ tC ∫ 1 0 |ζ 0 (t/ε)| dη ≤ tCe -M 0 t/ε (80) Iestimate
We used that d dt (-αφ(u ε (ηt)) + Φ(u ε (ηt), φ(u ε (ηt)))) = 0 and Assumption 4. Insertion of (80) into (78) yields

ε d dt V ε = -(α - ∂Φ ⋆ ∂v )V ε + g 1 (t) with g 1 (t) := -ε ∇ u Φ(u ε , φ(u ε )) α -∂Φ ∂v (u ε , φ(u ε )) • f (u ε , v ε ) + I(t), satisfying |g 1 | ≤ C(te -M 0 t/ε + ε) ≤ Cε.
Assumption 4, Corollary 1 and Remark 3 imply -(α -∂Φ ⋆ ∂v ) < -M 0 < 0. Using Young's inequality we obtain

|V ε (t)| ≤ ∫ t 0 e -M 0 (t-τ )/ε |g 1 (τ )| ε dτ ≤ Cε. ( 81
) Y1

Proof of Proposition 2

Sec5.2

We obtain

dδ ε dt = f (A + δ ε , φ(A + δ ε , t) + ζ 0 ( t ε ) + V ε , t) -f (A, φ(A, t), t) δ ε (0) = 0. ( 82
)
Due to Lipschitz continuity of the right hand side there exists a unique solution to this problem. It holds

f (A + δ ε , φ(A + δ ε , t), t) -f (A, φ(A, t), t) = (∇ u f + ∂ v f ⊗ ∇φ)(η(t))δ ε (83) |f (A + δ ε , φ(A + δ ε , t) + ζ 0 ( t ε ) + V ε , t) -f (A + δ ε , φ(A + δ ε , t), t)| ≤ C(e -M 0 t/ε + ε), ( 84
)
where η(t) is between A(t) and A(t) + δ ε (t). For the latter estimate we use Proposition 1 and (41). We obtain

dδ ε dt = (∇ u f + ∂ v f ⊗ ∇φ)(η(t))δ ε + g 2 (t)
with |g 2 | ≤ C(e -M 0 t/ε + ε). We note that if ε is small enough, η remains in a predefined vicinity of A. Therefore, we can use the exponential dichotomy property from Remark 2. Then Assumption 5 implies

|δ ε | ≤ Cε||e -M t || L 1 (R) ≤ Cε.
(85) deltaep1s

Proof of Remark 6

Sec5.2A

We start from equation

dδ ε dt = (∇ u f + ∂ v f ⊗ ∇φ)(η(t))δ ε + g 2 (t). ( 86
) eqdelta It implies d|δ ε | dt -(∇ u f + ∂ v f ⊗ ∇φ)(η(t))δ ε • δ ε |δ ε | ≤0 = g 2 (t) • sign δ ε . (87) estdelta1
The integration of (87) yields estimate (59).

Proof of Proposition 3

Sec5.3

We set ξ(x) := α -∂Φ ∂v (x, φ(x)) and

U ε = v ε -φ(u ε ) + ε ξ(u ε ) d dt φ(u ε ) -ζ 0 (t/ε) -εζ 1 (t/ε),
where ζ 0 is given by (40) and ζ 1 by (62).

Then, U ε (0) = 0 and |ζ 1 (t)| ≤ Ce -αt . Using (40) and definition of φ we obtain the following equation for

U ε ε d dt U ε = -αv ε + Φ(u ε , v ε ) -ε d dt φ(u ε ) + ε 2 d dt ( d dt φ(u ε ) ξ(u ε ) ) -ε d dt ζ 0 (t/ε) -ε 2 d dt ζ 1 (t/ε) = -αU ε -αφ(u ε ) -αζ 0 (t/ε) -αεζ 1 (t/ε) + ε ∂Φ ∂v (u ε , φ(u ε )) d dt φ(u ε ) ξ(u ε ) +Φ(u ε , v ε ) + ε 2 d dt ( d dt φ(u ε ) ξ(u ε ) ) -ε d dt ζ 0 (t/ε) -ε 2 d dt ζ 1 (t/ε) = -(α -∂ v Φ (1) )U ε + Φ ( u ε , φ(u ε ) - d dt φ(u ε ) ξ(u ε ) ε + ζ 0 (t/ε) + εζ 1 (t/ε) ) -Φ(u ε , φ(u ε )) + Φ(u 0 , φ(u 0 )) -Φ(u 0 , φ(u 0 ) + ζ 0 (t/ε)) +ε 2 d dt ( d dt φ(u ε ) ξ(u ε ) ) + ε ∂ v Φ(u ε , φ(u ε )) ξ(u ε ) d dt φ(u ε ) = -(α -∂ v Φ (1) )U ε + G ε + ε 2 d dt ( 1 ξ(u ε ) d dt φ(u ε ) ) . ( 88
) caluODE Here, ∂ v Φ (1) := ∂ v Φ(u ε , v (1) ) with v (1) between φ(u ε ) -ε ξ(uε) d dt φ(u ε ) + ζ 0 (t/ε) + εζ 1 (t/ε) and φ(u ε ) -ε ξ(uε) d dt φ(u ε ) + ζ 0 (t/ε) + εζ 1 (t/ε) + U ε . Next, we estimate the term G ε . Φ ( u ε , φ(u ε ) - ε ξ(u ε ) d dt φ(u ε ) + ζ 0 (t/ε) + εζ 1 (t/ε) ) -Φ(u ε , φ(u ε )) +Φ(u 0 , φ(u 0 )) -Φ(u 0 , φ(u 0 ) + ζ 0 (t/ε)) + ε ∂ v Φ(u ε , φ(u ε )) ξ(u ε ) d dt φ(u ε ) = Φ(u 0 , φ(u 0 )) -Φ(u 0 , φ(u 0 ) + ζ 0 (t/ε)) -Φ(u ε , φ(u ε )) +Φ ( u ε , φ(u ε ) + ζ 0 (t/ε) + εζ 1 (t/ε) ) + ∂ vv Φ (2) ε 2 2(ξ ( u ε )) 2 ( d dt φ(u ε ) ) - ( ∂ v Φ ( u ε , φ(u ε ) + ζ 0 (t/ε) + εζ 1 (t/ε) ) -∂ v Φ(u ε , φ(u ε )) ) ε ξ(u ε ) d dt φ(u ε ), (89) caluestim
where ∂ vv Φ (2) = ∂ vv Φ(u ε , v (2) ) for an appropriate value v (2) .

We start by estimating the two residuals in (89). Due to the boundedness of derivatives of Φ and f (u ε , v ε ), and exponential decay of the initial layers we obtain

( ∂ v Φ ( u ε , φ(u ε ) + ζ 0 (t/ε) + εζ 1 (t/ε) ) -∂ v Φ(u ε , φ(u ε )) ) d dt φ(u ε ) ξ(u ε ) ε + ∂ vv Φ (2) 2(ξ(u ε )) 2 ( d dt φ(u ε ) ) ε 2 ≤ C ( ε 2 + εe -κt/ε ) . ( 90
) caluestim
Estimating the remaining terms in (89) involves the boundary layers ζ 0 and ζ 1 and the use of the integral remainder for Taylor's formula:

|Φ ( u ε , φ(u ε ) + ζ 0 (t/ε) + εζ 1 (t/ε) ) -Φ(u ε , φ(u ε )) +Φ(u 0 , φ(u 0 )) -Φ(u 0 , φ(u 0 ) + ζ 0 (t/ε))| = ∫ 1 0 (u ε -u 0 ) • ∇ x Φ ( x, φ(x) + ζ 0 (t/ε) ) x=u 0 +η(uε-u 0 ) dη - ∫ 1 0 (u ε -u 0 ) • ∇ x Φ ( x, φ(x) ) x=u 0 +η(uε-u 0 ) dη + Cεe -κt/ε ≤ t m ∑ i=1 d dt u i ε ∞ ∂ v ∇ x i Φ ( x, φ(x) ) ∞ Ce -κt/ε + Cεe -κt/ε ≤ tCe -κt/ε + Cεe -κt/ε . ( 91
(A) = Φ(A, φ(A)) yield d dt v ε (t) = α ε -v ε (t) + φ(A(t)) - Φ(A(t), φ(A(t))) α + Φ(u ε (t), v ε (t)) α ≤ C(1 + 1 ε e -M 0 t/ε ). ( 92
) derivvest Using ∇φ(x) = ∇ u Φ(x, φ(x)) ξ(x)
, we obtain

g 4 (t) := d dt ( 1 ξ(u ε ) d dt φ(u ε ) ) = d dt ( ∇Φ(u ε , φ(u ε )) • f (u ε , v ε ) (ξ(u ε )) 2 ) = ∇ u ( ∇Φ(u ε , φ(u ε )) • f (u ε , v ε ) (ξ(u ε )) 2 ) • f (u ε , v ε ) + ∂ ∂v ( ∇Φ(u ε , φ(u ε )) • f (u ε , v ε ) (ξ(u ε )) 2 ) d dt v ε . ( 93
) caluestim
Due to the boundedness of derivatives of f and Φ and equation (92) it follows

|g 4 (t)| ≤ C ( 1 + 1 ε e -M 0 t/ε ) . ( 94
) caluestim
Inserting equations ( 90) and (91) in equation ( 89), and finally equation (89) and estimate (94) in equation (88), we obtain

ε d dt U ε = -(α -∂ v Φ (1) )U ε + g 5 (t) (95) with |g 5 (t)| ≤ C(ε 2 + (ε + t)e -κt/ε ).
We obtain

d dt δε = f ( A + εu 1 + εψ 1 (t/ε) + δε , U ε +φ(A(t) + εu 1 (t) + εψ 1 (t/ε) + δε ) -ε ∇φ(u ε ) • f (u ε , v ε ) ξ(u ε ) +ε ∇φ(A) • d dt A ξ(A) -ε ∇φ(A) • d dt A ξ(A) + ζ 0 (t/ε) + εζ 1 (t/ε) ) - d dt u 1 ε -ε d dt ψ 1 (t/ε) -f (A, φ(A)) = f ( A + εu 1 , -ε ∇φ(A) • d dt A ξ(A) + ζ 0 (t/ε) + φ(A + εu 1 ) ) +D u f * (εψ 1 (t/ε) + δε ) - d dt u 1 ε -ε d dt ψ 1 (t/ε) -f (A, φ(A)) +∂ v f * ( U ε -ε ∇φ(u ε ) • f (u ε , v ε ) ξ(u ε ) + ε ∇φ(A) • d dt A ξ(A) + εζ 1 (t/ε) ) , ( 103
)
deltau2OD
where ∂ v f * = ∂ v f (u ε , v * ) for an appropriate value v * and

D u f * = ∇ u f ( x, φ(x) - ε ξ(A) ∇φ(A) • d dt A + ζ 0 (t/ε) ) x=u * +∂ v f ( x, φ(x) - ε ξ(A) ∇φ(A) • d dt A + ζ 0 (t/ε) ) ⊗ ∇φ(x) x=u *
for an appropriate value u * between u ε and A + εu 1 .

The perturbation ε ξ(A)

∇φ(A)• d dt A+ζ 0 (t/ε) tends to zero for ε tending to zero and t ≥ O [START_REF] Agarwal | An Introduction to Ordinary Differential Equations[END_REF].

Therefore, D u f * satisfies the exponential dichotomy from Assumption 5. Now we estimate the last term in (103): )

U ε -ε ∇φ(u ε ) • f (u ε , v ε ) ξ(u ε ) + ε ∇φ(A) • d dt A ξ(A) = U ε -ε ∇φ(u ε ) • f (u ε , v ε ) ξ(u ε ) + ( ∇φ(A) • f (A, φ(A)) ξ(A) - ∇φ(u ε ) • f (u ε , φ(u ε )) ξ(u ε ) ) ε + ∇φ(u ε ) • f (u ε , φ(u ε )) ξ(u ε ) ε = Ûε + ( ∇φ(A) • f (A, φ(A)) ξ(A) - ∇φ(u ε ) • f (u ε , φ(u ε )) ξ(u ε ) ) ε (104) remainder with Ûε := U ε -ε ∇φ(u ε ) ξ(u ε ) • (f (u ε , v ε ) -f (u ε , φ(u ε ))). ( 105
+ D u f * δε - d dt u 1 ε -ε d dt ψ 1 (t/ε) -f (A, φ(A)) + ∂ v f * Ûε,2 + ε(G 1 + G 2 ). ( 110
) -f (A, φ(A)) + D u f * δε + [ D u f A u 1 -∂ v f A ∇φ(A) • d dt A ξ(A) ε - d dt u 1 ] ε -ε d dt ψ 1 (t/ε) + ∂ v f * Ûε,2 + ε(G 5 + G 6 ) (

  ) DecayEta where C and κ are positive constants. Consequently v 0 (τ ) = φ(A) + O(e -κτ ).

Figure 1 :

 1 Figure 1: The tubular domain T b,λ .

(A

  ) delatu2ODUsing Assumption 5 we obtain the estimates||G 1 || L ∞ (R + ) ≤ Cε and |G 2 (t)| ≤ Ce -M 0 t/ε . Using boundedness of u 1 , we calculate f + εu 1 , -ε ∇φ(A) • d dt A ξ(A) + ζ 0 (t/ε) + φ(A + εu 1 ) 0 (t/ε) + φ(A) ) + D u f A εu 1 -∂ v f A ∇φ(A) • d dt A ξ(A) ε +ε(G 3 + G 4 ) (111) fApprox with ||G 3 || L ∞ (R + ) ≤ Cε and |G 4 (t)| ≤ Ce -M 0 t/ε . Insertion of equation (111) into equation (110) for δε leads to d dt δε = f ( A, ζ 0 (t/ε) + φ(A)

1 0 (

 10 112) deltau2OD with ||G 5 || L ∞ (R + ) ≤ Cε and |G 6 (t)| ≤ Ce -M 0 t/ε .Finally, insertion of equations (101) and (100) into (112) leads tod dt δε = f ( A, ζ 0 (t/ε) + φ(A) ) -f (A, φ(A)) + D u f * δε -f (u 0 , φ(u 0 ) + ζ 0 (t/ε)) + f (u 0 , φ(u 0 )) + ε(G 5 + G 6 ) + ∂ v f * Ûε,2 . (113) deltau2ODA precise Taylor estimate results in||f (A, ζ 0 (t/ε) + φ(A)) -f (A, φ(A)) -f (u 0 , φ(u 0 ) + ζ 0 (t/ε)) + f (u 0 , φ(u 0 ))|| ∞ = ∫ ∇ u f (u, φ(u) + ζ 0 (t/ε)) | u=u 0 +ν(A-u 0 ) -D u f (u, φ(u)) | u=u 0 +ν(A-u 0 ) ) (A -u 0 )dν ∞ ≤ C max{ε 2 , te -M 0 t/ε } d dt A ∞ ∥∂ v D u f (u, φ(u) + v)∥ ∞ =: G 7 max{ε 2 , te -M 0 t/ε }. (114) estIn the final step we use estimate (114) and the exponential dichotomy for D u , to obtain the inequality| δε | ≤ ∫ t 0 e -M (t-τ ) (|εG 5 (τ ) + εG 6 (τ ) + ∂ v f * Ûε,2 | + τ e -M 0 τ /ε G 7 ) dτ ≤ Cε 2 + C|| Ûε,2 || L 1 (R + ) ≤ Cε 2 . (115)

  ) caluestim

	It remains to estimate the term ε 2 d dt	(	1 ξ(uε)	) dt φ(u ε ) d	in (88). We start by the estimate for

d dt v ε . Theorem 1 and the equation αφ

  ) calutilde Since | Ûε | ≤U ε + εC|V ε + ζ 0 (t/ε)|, || L ∞ (R + ) ≤ Cε 2 , || Ûε,2 || L 1 (R + ) ≤ Cε 2 .

								(106)
	using a similar argument as in (97), we obtain
					Ûε =: Ûε,1 + Ûε,2	(107) calutilde
	with
				|| Ûε,1 (108)
	Next, Theorem 1 implies
		(	∇φ(A) • f (A, φ(A)) ξ(A)	-	∇φ(u ε ) • f (u ε , φ(u ε )) ξ(u ε )	)	≤ Cε.	(109) remainder
	Insertion of (109) and (107) into (104) and insertion of (104) into equation (103) yield
	d dt	δε = f	( A + εu 1 , -ε	∇φ(A) • d dt A ξ(A)	+ ζ

0 (t/ε) + φ(A + εu 1 )

Acknowledgements

AM-C and TS are supported by the Collaborative Research Center, SFB 873 "Maintenance and Differentiation of Stem Cells in Development and Disease" from the German Research Council (DFG). Research visits of AM at the Heidelberg University were supported by the Romberg Professorship of the Heidelberg Graduate School, Interdisciplinary Center for Scientific Computing (IWR), Heidelberg University.

th2 Theorem 2 Let Assumptions 2-5 be satisfied. For ε small enough it holds

where A RG fulfills the 1st order RG equation

Proof. First we observe that

Since εψ 1 ( t ε ) decays exponentially, it is enough to estimate the error function d ε : = A RG -A -εu 1 . Using the ODEs for A RG , A and εu 1 , we obtain

where G * is bounded. Error estimate (i) now follows directly.

In order to prove (ii), we use Proposition 3 and the estimate (i).

Example of an application in biology

Sec3

We consider the following system of ordinary differential equations

where p, k, d are positive constants, a ∈ (0.5, 1) and u 0 1 , u 0 2 , v 0 are positive. This system describes time evolution of a cell system maintained by a stem cell population that differentiates to mature cells, e.g., blood cells [START_REF] Marciniak-Czochra | Modeling asymmetric cell division in hematopietic stem cells regulation of self-renewal is essential for efficient repopulation[END_REF], [START_REF] Stiehl | Characterization of stem cells using mathematical models of multistage cell lineages[END_REF]. Then v ε is interpreted as the concentration of signaling molecules and u ε,1 , u ε,2 as counts of stem and committed cells. The existence of (locally) stable positive equilibria and the dynamics in the neighborhood of such equilibria are of biological interest [START_REF] Nakata | Stability analysis of multicompartment models for cell production systems[END_REF], [START_REF] Stiehl | Characterization of stem cells using mathematical models of multistage cell lineages[END_REF]. The 0th-order RG approximation ("quasi steady state approximation") of this model Using the explicit formula for U ε and Young's inequality yields

where

(97) mathcalUs

Proof of Proposition 4

Sec5.4

As above, we set ξ(x) := α -∂Φ ∂v (x, φ(x)). Adopting the notations from Proposition 3, we obtain

In order to obtain an O(ε 2 ) approximation, we expand u ε as

Then it is required that u 1 (0) + ψ 1 (0) = 0. We obtain Note that according to Assumptions 3 and 5 u 1 is bounded. We now estimate the approximation error δε (t) = u ε (t) -A(t) -εu 1 (t) -εψ 1 (t/ε).

(102) deltaudef