MD, PhD Alice Courties

MD Ariane Do

Sarah Leite

Manon Legris

PhD Laure Sudre

Audrey Pigenet

MD Juliette Petit

MD, PhD Geoffroy Nourissat

MD Adeline Cambon-Binder

PhD Uwe Maskos

MD, PhD Francis Berenbaum
email: francis.berenbaum@aphp.fr

MD, PhD Jérémie Sellam

Prof Berenbaum Francis

The role of the non-neuronal cholinergic system in inflammation and degradation processes in osteoarthritis Running Head : Non-neuronal cholinergic system involvement in osteoarthritis

Objectives: The non-neuronal cholinergic system represents non-neuronal cells that have the biochemical machinery to synthetize de novo and/or respond to acetylcholine (ACh). We investigated this biochemical machinery in chondrocytes and its involvement in osteoarthritis (OA).

Methods: Expression of the biochemical machinery for ACh metabolism and nicotinic ACh receptors (nAChR), particularly α7nAChR, in human OA and murine chondrocytes was determined by PCR and ligand-binding. We investigated the mRNA expression of the human duplicate α7nACh subunit, called CHRFAM7A, which is responsible for a truncated α7nAChR.

We assessed the effect of nAChR on chondrocytes activated by IL1β and the involvement of α7nAChR using chondrocytes from wild-type (WT) and α7-deficient Chrna7 -/-mice. The role of α7nAChR in OA was explored after medial meniscectomy (MNX) in WT and Chrna7 -/-mice.

Results: Human and murine chondrocytes express the biochemical partners of the non-neuronal cholinergic system (n=5) and a functional α7nAChR at their cell surface. The expression of CHRFAM7A in human OA chondrocytes correlated positively with MMP3 (r=0.38, p<0.05) and MMP13 (r=0.48, p<0.05) expression (n=23). Nicotine decreased the IL1β-induced IL6 and MMPs expression in a dose-dependent manner in WT but not in Chrna7 -/-chondrocytes. MNX Chrna7 -/- mice displayed more severe OA cartilage damage (OARSI score of 4.46±1.09, n=7) than MNX WT mice (3.05±0.9, n=9, p<0.05).

INTRODUCTION

While osteoarthritis (OA) is the most common musculoskeletal disease, available treatments still only treat symptoms and are either not highly effective or associated with safety issues. Therefore, it remains an unmet need to develop new therapeutic targets based on a better understanding of OA pathophysiology.

The interplay between the nervous system and inflammation has been mainly described during infection or acute inflammation [START_REF] Tracey | Neurons Are the Inflammatory Problem[END_REF]. Along this line, the autonomic nervous system, particularly the parasympathetic system, has shown immunomodulatory and anti-inflammatory properties after stimulation of the vagus nerve [START_REF] Borovikova | Vagus nerve stimulation attenuates the systemic inflammatory response to endotoxin[END_REF]. The parasympathetic system is composed of cholinergic neurons that release acetylcholine (ACh) as the chemical mediator of neurotransmission. However, ACh is not only a neurotransmitter, as it can be synthesized by non-neuronal cells such as T cells or resident cells (e.g., keratinocytes) [START_REF] Wessler | Acetylcholine beyond neurons: the non-neuronal cholinergic system in humans[END_REF]. These cells can also express the biochemical machinery responsible for the transport and degradation of synthesized ACh and nicotinic and/or muscarinic ACh receptors. This ACh production and responsive system (i.e., ACh receptors) define the non-neuronal cholinergic system.

Anti-inflammatory properties of non-neuronal ACh produced by T cells have been demonstrated

in sepsis, viral infections and rheumatoid arthritis involving activation of the alpha 7 nicotinic ACh receptor (α7nAChR) expressed by splenic macrophages [START_REF] Borovikova | Vagus nerve stimulation attenuates the systemic inflammatory response to endotoxin[END_REF][START_REF] Cox | Choline acetyltransferase-expressing T cells are required to control chronic viral infection[END_REF][START_REF] Wang | Nicotinic acetylcholine receptor alpha7 subunit is an essential regulator of inflammation[END_REF]. The α7nAChR is composed of five α7 subunits (Chrna7). The α7nAChR mediates its biological effect through ionic flux, especially intracellular calcium entrance. However, it has been suggested that ACh could activate direct metabotropic signalling, such as JAK2-STAT3 or MAP-kinase signalling [START_REF] Marrero | Convergence of alpha 7 nicotinic acetylcholine receptor-activated pathways for anti-apoptosis and anti-inflammation: central role for JAK2 activation of STAT3 and NF-kappaB[END_REF]. Recently in humans, genetic rearrangement between a duplicate of the 5-10 exons of Chrna7 and a copy of the FAM7A gene led to a specific human duplicate responsible for the production of a truncated protein (CHRFAM7A) [START_REF] Gault | Genomic organization and partial duplication of the human alpha7 neuronal nicotinic acetylcholine receptor gene Accepted Article (CHRNA7)[END_REF]. This truncated protein may decrease the biological response to ACh and thus promote inflammation [START_REF] Costantini | A Human-Specific α7-Nicotinic Acetylcholine Receptor Gene in Human Leukocytes: Identification, Regulation and the Consequences of CHRFAM7A Expression[END_REF][START_REF] Costantini | Uniquely human CHRFAM7A gene increases the hematopoietic stem cell reservoir in mice and amplifies their inflammatory response[END_REF].

Accepted Article

This article is protected by copyright. All rights reserved During OA, chondrocytes are activated under different stresses (mechanical, inflammatory, metabolic) and are then responsible for the production of inflammatory and catabolic mediators that induce cartilage degradation. Two recent publications showed that systemic nicotinic ACh receptor activation by intraperitoneal injection of nicotine may limit OA lesions in a monoiodoacetate murine model of OA [START_REF] Teng | Nicotine Attenuates Osteoarthritis Pain and Matrix Metalloproteinase-9 Expression via the α7 Nicotinic Acetylcholine Receptor[END_REF][START_REF] Liu | Activation of α7 nicotinic acetylcholine receptors prevents monosodium iodoacetate-induced osteoarthritis in rats[END_REF]. Using a chemical antagonist compound (i.e., methyllycaconitine), they suggested that α7nAChR was involved in macrophage and chondrocyte activation. However, these results were not confirmed in another study [START_REF] Bock | What is the effect of nicotinic acetylcholine receptor stimulation on osteoarthritis in a rodent animal model?[END_REF].

Likewise, the role of nicotinic ACh receptors in OA remains controversial. Furthermore, it is unknown whether articular chondrocytes display the non-cholinergic system to produce ACh and express nicotinic ACh receptors or whether activation of cholinergic receptors influences the chondrocytic production of cytokines and metalloproteases. The objective of this study was to determine the expression and role of the non-neuronal cholinergic system on chondrocyte activation and OA, especially through the involvement of α7nAChR.

MATERIALS AND METHODS

Animals. α7-deficient Chrna7 -/-(B6.129S7-Chrna7 tm1Bay /J) and wild-type (WT) C57BL/6 mice were used for in vitro and in vivo experiments [START_REF] Orr-Urtreger | Mice deficient in the alpha7 neuronal nicotinic acetylcholine receptor lack alpha-bungarotoxin binding sites and hippocampal fast nicotinic currents[END_REF]. Animals were fed ad libitum. Induction of experimental OA. Meniscectomy (MNX) was performed on 10-12-week-old male Chrna7 -/-and WT mice as previously described [START_REF] Kamekura | Osteoarthritis development in novel experimental mouse models induced by knee joint instability[END_REF][START_REF] Kadri | Osteoprotegerin inhibits cartilage degradation through an effect on trabecular bone in murine experimental osteoarthritis[END_REF]. Briefly, MNX was performed by a trained surgeon : after cutaneous incision, the medial meniscotibial ligament was sectioned and medial meniscus extrusion and sectioning was performed on the right knee. Mice were euthanized 8 weeks after MNX for analyzes. Sham procedure was performed on the right knee of independent male WT and Chrna7 -/-mice through the same approach but without any ligament transection or meniscectomy.

Histological analysis. After sacrifice, the right knee joints were removed and fixed for 24 hours in 4% paraformaldehyde-phosphate-buffered saline (PBS) (PFA) solution at 4°C. Knee joints were

Accepted Article

This article is protected by copyright. All rights reserved then decalcified in 500 mM EDTA for 4 weeks, and all tissues were embedded in paraffin. Sagittal sections of 5 µm were obtained. Knee tissue sections were stained with Safranin O/fast green for OARSI histologic OA scoring at the tibial and femoral condyles [START_REF] Glasson | The OARSI histopathology initiative -recommendations for histological assessments of osteoarthritis in the mouse[END_REF]. For each knee, the final score is the mean score of tibial and femoral condyles obtained from 3 slides separated by 50 µM by two readers blinded to the genotype. Synovitis scores obtained were the sum of the score of the synovial hyperplasia (0-3) and of sub-synovial inflammation (0-3) as previously published [START_REF] Jackson | Depletion of protease-activated receptor 2 but not protease-activated receptor 1 may confer protection against osteoarthritis in mice through extracartilaginous mechanisms[END_REF].

Osteophytes scoring was performed as previously published on a 0-3 scale [START_REF] Kamekura | Osteoarthritis development in novel experimental mouse models induced by knee joint instability[END_REF].

Alpha-bungarotoxin staining. Human and WT murine chondrocytes were grown on slides for 5 days and incubated with biotinylated-bungarotoxin (Thermo Fisher Scientific, Waltham, Massachusetts, USA) in the cell culture medium at 4°C overnight at 5 µg/L and 20 µg/L, respectively. Cells were washed three times with PBS and then fixed for 15 min at room temperature in 4% PFA. Slides were then incubated at room temperature in streptavidinhorseradish peroxidase (HRP) (Vectastain ® -Vector Laboratories, Inc. CA, USA) solution using DAB (Vector Laboratories, Inc. CA, USA) substrate and counterstained with haematoxylin. Semiquantitative scoring was performed on four different zones of 0.1 mm 2 of the slides, and the stained cells were reported relative to the total number of cells.

Primary culture of murine articular chondrocytes. Murine chondrocytes were obtained from the femoral head and the knee articular cartilage of 5-to 6-day-old Chrna7 -/-and WT mice, as previously described [START_REF] Gosset | Primary culture and phenotyping of murine chondrocytes[END_REF]. Briefly, after matrix digestion, chondrocytes were seeded at 8 × 10 3 cells per cm 2 and expanded for a week in Dulbecco's modified Eagle's medium (DMEM) with 1 g/L glucose supplemented with 100 U/mL penicillin (P), 100 mg/mL streptomycin (S), 4 mM glutamine (Glu) and 10% foetal calf serum (FCS). Each littermate was used for one experiment. Cells were fasted for 24 hours before treatment in serum-free media DMEM P/S/Glu containing 0.1% bovine serum albumin (BSA).

Treatment of primary cultures of murine chondrocytes. Chrna7 -/-and WT chondrocytes were stimulated for 24 h with IL-1β (10 ng/mL) (PeproTech, Rocky Hill, NJ, USA) with or without pretreatment for 15 minutes with nicotine (Sigma, Saint-Louis, USA) at 1, 10 or 100 µM. After 24 h, cell lysates were collected for mRNA isolation, and media were collected for protein assessment by enzyme-linked immunosorbent assay (ELISA).

Accepted Article

This article is protected by copyright. [START_REF] Grønlien | Distinct profiles of alpha7 nAChR positive allosteric modulation revealed by structurally diverse chemotypes[END_REF]. After 24 h, cell lysates were collected for mRNA analysis.

Calcium assay. Primary cultures of murine chondrocytes in 96-well plates were starved after 1 week of culture. Measurements were performed using a Fluo-4 NW calcium assay kit (Molecular Probes, Eugene, USA) following the manufacturer's instructions. Briefly, after being rinsed with saline phosphate buffer (PBS), the cells were incubated with probenecid in the kit buffer to block cationic organic transporters for 30 minutes at 37° and then 30 minutes at room temperature. Then the live analysis until 3.8 minutes was performed using an inverted fluorescence an inverted microscope (Olympus IX83 inverted), immediately after the addition of an agonist: ionomycin at 10 M (positive control) or nicotine at 10 M. The fluorescence was then quantified by Image J.

RNA extraction and RT-PCR.

Total RNA was extracted from chondrocytes using an RNeasy kit (Qiagen, Hilden, Germany), and concentrations were determined using a Nanodrop (DeNovix DS-11). Reverse transcription involved 500 ng total RNA with an Omniscript RT kit (Qiagen, Hilden, Germany). mRNA levels of interleukin 6 (IL6), metalloproteases (MMP)-3, MMP9, MMP13, nicotinic ACh receptor subunits (α1-7,9 and β 1-4) and the human duplicate CHRFAM7A were quantified using a Light Cycler LC480 (Roche Diagnostics) by qRT-PCR. Levels of mRNA were normalized to those of murine hypoxanthine guanine phosphoribosyltransferase (HPRT) and to human 18S. The relative quantitative value is equal to 2 -CT where CT is equal to (CT target-CT HPRT or 18S). For nAchR expression, we therefore normalized on a calibrator (a mix of mRNA samples).

Accepted Article

RT-PCR was performed to determine the expression of factors in the cholinergic system, i.e., carnityl acetyltransferase (carAT), choline acetyltransferase (ChAT), vesicular acetylcholine transporter (VAChT), acetylcholine esterase (AChE), butyrylcholine esterase (BChE) and choline transporter-1 (CHT1), in human and murine chondrocytes. For RT-PCR, PCR was performed with a final volume of 20 µL containing 10 µL of Taq polymerase (Qiagen), 0.25 µM of each forward and reverse primer (Supplementary S1 table 1 for murine and table 2 for human primers) and 6 µL of cDNA. For each primer, we determined the best melting temperature and applied a negative control (water) and a positive control (murine and human cDNA from the spinal cord).

Amplification was performed for 40 cycles. For electrophoretic separation, RT-PCR products were subjected to 2 % TRIS-acetate-EDTA agarose gel electrophoresis.

ELISA. The protein concentrations of IL6 and MMP3 released in the supernatant by WT and

Chrna7 -/-chondrocytes after IL1β with or without nicotine were measured using a Quantikine ELISA kit (R&D Systems, Lille, France).

Statistical analysis.

A non-parametric Wilcoxon test and a Mann-Whitney test were used for matched paired and unpaired data, respectively. A P-value <0.05 indicated statistical significance.

Data are expressed as the mean±standard error of the mean (SEM). Statistical analyzes were performed using GraphPad Prism 5 (GraphPad Software, San Diego, CA, USA).

RESULTS

Murine and human OA chondrocytes express the non-neuronal cholinergic system

Human OA and WT murine chondrocytes expressed carAT, VAChT, AChE (mRNA) in 5 independent experiments using 5 different samples. All murine chondrocytes expressed CHT1 while none of the human chondrocytes expressed it. Four murine chondrocytes and two human chondrocytes expressed BChE, and none of them expressed ChAT. The mRNA expression of each molecular actor of the non-neuronal cholinergic system was positive for cDNA in the spinal cord (i.e., positive control) and negative in water, confirming the specificity of the primers (Figure 1A).

Accepted Article

Activation of nAChR exerts anti-inflammatory and anti-catabolic action on murine and human chondrocytes

We focused on the role of nicotine, which binds all the so-called nicotinic receptors, on chondrocyte activation. We incubated murine WT chondrocytes with nicotine stimulated by IL1β.

Nicotine decreased the mRNA expression of IL6 in a dose-dependent manner between 1 and 10 M by 37.3±8%; 70.8±6% and 46.8±10% for 1, 10 and 100 µM nicotine, respectively (n=5, p<0.05) (Figure 1B). Additionally, nicotine decreased the mRNA expression of MMP3 in a dosedependent manner between 1 and 10 M (by 46.5±8%; 68.3±10% and 67.8±11% for 1, 10 and 100 µM nicotine, respectively, p<0.05) (Figure 1B) and of MMP13 and MMP9 mRNA (Supplementary S2). The anti-inflammatory and anti-catabolic effects reached a plateau between the 10 and the 100 µM dose. These results were confirmed at the protein level since nicotine decreased the levels of IL1β-induced IL6 in a dose-dependent manner according to ELISA with decreases of 38.9±16.9% (not significant), 65.6±4%, and 67.7±6.4% for 1, 10 and 100 µM nicotine, respectively (n=5; p<0.05) (Figure 1C). Nicotine also decreased the levels of MMP3 secreted by murine chondrocytes after IL1β stimulation, with a decrease of 54±11% for 100 µM nicotine (n=5; p<0.05) (Figure 1C).

Murine and human chondrocytes express several nicotinic receptor subunits, especially the functional α7nAchR

Since we identified anti-inflammatory and anti-catabolic effects of nicotine on chondrocytes, we questioned which nAchR subunits are expressed by human and murine chondrocytes. Using qPCR, we found that murine chondrocytes expressed α4, α5, α6, α7 and β4 subunits (between 31 and 33 cycles), while human OA chondrocytes expressed α5, α6, α7 and β2, β4 subunits (between 30 and 35 cycles) (Figure 2A). As both chondrocytes expressed the Chrna7 subunit, we determined whether murine and human chondrocytes express the functional homopentameric receptor α7nAchR using α-bungarotoxin, a specific antagonist of functional α7nAChR. Murine and human chondrocytes bound α-bungarotoxin in two independent experiments each (Figure 2B and 2C). Using semi-quantitative scoring, we found that 27.8±4.8% of murine chondrocytes and 36.3±4.5% of human OA chondrocytes expressed a functional homopentamer α7nAchR (Figure 2B and2C).

Accepted Article

This article is protected by copyright. All rights reserved

Human chondrocytes express both native Chrna7 and duplicate CHRFAM7A, and

CHRFAM7A expression positively correlates with MMP3 and -13 expression

To better decipher α7nAChR expression in human cells, we evaluated the expression of both native Chrna7 and the human-specific duplicate CHRFAM7A subunit in human OA chondrocytes.

The structural differences between native Chrna7 and duplicate CHRFAM7A are explained in Figure 3A. We evaluated the expression of both forms in human OA chondrocytes from 23 patients with knee OA undergoing arthroplasty (69% women, mean ± standard deviation 69.6±7.9 years old). Human OA chondrocytes expressed both the duplicate CHRFAM7A and the native Chrna7 subunits, but the mRNA expression of CHRFAM7A was significantly higher than that of native Chrna7 (p=0.046) (Figure 3B). Interestingly, basal human OA chondrocyte CHRFAM7A mRNA expression correlated with MMP3 (r=0.38, n=23, p<0.05) and MMP13 (r=0.48, n=23, p<0.05) mRNA expression (Figure 3C). High-level expression of CHRFAM7A was thus associated with a more pronounced catabolic profile since the expression of MMP3 and MMP13 was significantly higher in chondrocytes expressing high levels of CHRFAM7A (defined by

CHRFAM7A≥median which was 1.4) compared with those who expressed low levels of CHRMAF7A (CHRFAM7A<median) (Figure 3D).

The α7nAChR pathway is critical for anti-inflammatory and anti-catabolic activity induced by cholinergic system activation

We hypothesized that the anti-inflammatory effect of nAChR activation was due to α7nAChR.

Therefore, we investigated the effect of nicotine on IL1-stimulated Chrna7 -/-chondrocytes. In contrast to the result in WT chondrocytes, nicotine had no protective effect and thus did not decrease IL6 and MMP3 mRNA expression at any dose in Chrna7 -/-chondrocytes (Figure 4A).

This phenomenon was confirmed at the protein level (Figure 4B). Therefore, the antiinflammatory and anti-catabolic effects of nicotine depend on α7nAchR. We excluded any compensation effect of other subunits in the Chrna7 -/-chondrocytes because the expression level of the subunits of nAChR (α1-7,9 and β1-4) did not differ between WT and Chrna7 -/-chondrocytes (not shown). To determine if the effect of nicotine was calcium-dependent, we performed calcium entrance experiments in WT and Chrna7 -/-murine chondrocytes and found no intracellular entrance of calcium in either type of chondrocyte, meaning that activation of α7nAchR involves other pathways (Figure 4C-E). To confirm the anti-inflammatory and anti-catabolic effects of the α7nAChR pathway in human chondrocytes, we tested the effect of a specific α7nAChR agonist,

Accepted Article

This article is protected by copyright. All rights reserved PNU-282987 (100 µM), on IL1-stimulated chondrocytes: PNU-282987 decreased the mRNA expression levels of IL6 (-23.1±5.9%), MMP3 (-16.5±6.5%), and MMP13 (-39.4±4.1%) induced by IL1 alone (Figure 4F).

Impairment of the ACh/Chrna7 axis by deletion of Chrna7 aggravates cartilage degradation in murine OA

Finally, we determined the involvement of α7nAChR in vivo using the MNX model of OA in WT and Chrna7 -/-mice. We did not find any macroscopic phenotypic or behavioural differences between WT and Chrna7 -/-mice before sacrifice. The weight was similar between the groups (27.4±0.7g for WT sham mice, n=7; 28.3±0.5 for Chrna7 -/-sham mice, n=5; 29.2±0.6g for WT MNX mice, n=9 and 28.3±0.8g for Chrna7 -/-MNX mice, n=7; not significant) (Figure 5B). MNX induces significant lesions of OA in both WT and Chrna7 -/-mice compared to sham knees with more cartilage lesions, more osteophyte and more synovitis (p<0.05) (Figure 5A and5C). In the sham experiments, there were no difference in the OARSI score between WT and Chrna7 -/- (0.42±0.1 for WT sham mice versus 0.84±0.51 for Chrna7 -/-sham mice, p=0.11) as well as in the osteophyte and synovitis scores (Figure 5D and5E).

Eight weeks after MNX, Chrna7 -/-mice display more severe OA cartilage damage than WT mice as assessed by the OARSI score (3.05±0.9 for WT MNX versus 4.46±1.09 for Chrna7 -/-MNX on a 0-6 score, 6 being the worst, p<0.05) (Figure 5A and5C). However, we did not find any difference in synovitis and osteophytes scores between WT and Chrna7 -/-mice after MNX (Figure 5D and5E).

DISCUSSION

Our understanding of the role of the non-neuronal cholinergic system as a critical regulator of inflammation and immunity is growing but remains controversial and poorly investigated in OA pathophysiology [START_REF] Fujii | Physiological functions of the cholinergic system in immune cells[END_REF][START_REF] Courties | Role of the autonomic nervous system in osteoarthritis[END_REF]. Here, we showed that human and murine chondrocytes are part of the non-neuronal cholinergic system allowing local auto/paracrine action of ACh on joint cells.

Stimulation of the functionally expressed homopentamer α7nAChR decreased the inflammatory and catabolic activation of human and murine chondrocytes induced by IL1. Furthermore, α7nAChR has a protective role in OA since the lack of α7nAChR was associated with more severe cartilage degradation in the murine MNX model.

Accepted Article

This article is protected by copyright. All rights reserved

The role of the non-neuronal cholinergic system has been previously investigated in other cell types. Likewise, non-neuronal production of ACh by T-cells is essential to control viral infections (4). Additionally, the stimulation of ACh synthesis by T-cells could counteract the chronic and systemic inflammation in rheumatoid arthritis and in sepsis . However, very little is known about the non-neuronal cholinergic system in joint tissues, especially in cartilage. The synovial membrane expresses ChAT, the main enzyme responsible for the production of ACh and carAT, another enzyme for ACh synthesis [START_REF] Schubert | Expression of the non-neuronal cholinergic system in human knee synovial tissue from patients with rheumatoid arthritis and osteoarthritis[END_REF]. In addition, the synovium expresses AChE and BChE, which are responsible for the degradation of ACh and several nicotinic ACh receptors, such as α2-7 and 9 subunits (but not β nicotinic subunits) and muscarinic receptors M1, M3-M5 [START_REF] Schubert | Expression of the non-neuronal cholinergic system in human knee synovial tissue from patients with rheumatoid arthritis and osteoarthritis[END_REF][START_REF] Grimsholm | Unexpected finding of a marked non-neuronal cholinergic system in human knee joint synovial tissue[END_REF][START_REF] Beckmann | Expression of choline and acetylcholine transporters in synovial tissue and cartilage of patients with rheumatoid arthritis and osteoarthritis[END_REF].

Murine osteoblasts also express the entire spectrum of biochemical enzymatic machinery and most nicotinic receptors [START_REF] Sato | Functional role of acetylcholine and the expression of cholinergic receptors and components in osteoblasts. FEBS Accepted Article This article is protected by copyright[END_REF]. Here, we showed for the first time that human OA and murine chondrocytes express the entire cholinergic system relevant for ACh synthesis (carAT), transport (VAChT), and degradation (AChE). Chondrocytes express carAT and no ChAT which is concordant with previous study that found carAT is the main enzyme responsible for ACh synthesis in non-neuronal cells [START_REF] Tucek | The synthesis of acetylcholine in skeletal muscles of the rat[END_REF].We found that choline transporter-1 (CHT1) is expressed by murine but not by human OA chondrocytes. However, Beckmann et al. established that other choline transporters were expressed in cartilage (choline transporter-like CTL-1 to 5, the organic cation transporter (OCT) 1 and 3 and the OCT novel-3) (28). These transporters allow choline recapture by the cells after ACh degradation which may constitute a rate-limiting step for ACh synthesis.

ACh acts on both muscarinic and nicotinic receptors, which are widely expressed by non-neuronal cells, making the cholinergic system complex. Since nicotinic receptors are also involved in inflammation and tissue remodelling in non-neuronal cells, we investigated their expression in chondrocytes. Nicotinic receptors are composed of 5 alpha and/or beta subunits, which comprise a heterogeneous family of hetero-or homopentamers. Chondrocytes express several alpha and beta nicotinic subunits, and most are preserved between mice and humans. Indeed, WT murine chondrocytes express α4, α5, α6, α7, and β4 subunits, and human OA chondrocytes express α5, α6, α7, β2, and β4 subunits. Based on previously published data on other cell types, different combinations of nAChR pentamers can be expressed by chondrocytes [START_REF] Albuquerque | Mammalian nicotinic acetylcholine receptors: from structure to function[END_REF]. But importantly,

Accepted Article

This article is protected by copyright. All rights reserved human and murine chondrocytes both express the α7nAChR homopentamer, as confirmed by αbungarotoxin binding, which is widely used for this purpose [START_REF] Chen | The alpha-bungarotoxin-binding nicotinic acetylcholine receptor from rat brain contains only the alpha7 subunit[END_REF].

Our study demonstrated a regulatory effect of α7nAChR in chondrocyte activation induced by IL1. In vitro, activation of nicotinic receptors by nicotine limits the induction of inflammatory and catabolic factors by IL1. Our results are in accordance with those reported by Liu et al., who showed that nicotine applied to IL1-stimulated rat chondrocytes inhibited signalling pathways, mitogen-activated protein kinases (MAPKs) and nuclear factor-kappa B (Nf-B) proinflammatory pathways. However, they did not investigate the downstream effect of nicotinic receptor activation on cytokines and MMPs or which nicotinic receptors are involved in these processes. Here, we demonstrated that the anti-inflammatory and anti-catabolic effects of nicotine are due to α7nAChR activation since nicotine did not counteract IL1 stimulation in Chrna7 -/- chondrocytes. In contrast to Liu et al., we used Chrna7 -/-cells, which are more robust and more specific than chemical antagonist compounds such as methyllycaconitine (MLA), an α7nAChR antagonist. We did not find any action of nicotine on the entry of calcium, contrary to what has been previously published on human growth plate chondrocytes [START_REF] Kawakita | Nicotine acts on growth plate chondrocytes to delay skeletal growth through the alpha7 neuronal nicotinic acetylcholine receptor[END_REF]. The state of differentiation could explain this discrepancy. Despite calcium entry, α7nAChR could have a direct metabotropic effect by a direct interaction with some substrates of signalling pathways (Nf-B and MAPKs), as previously reported in chondrocytes [START_REF] Liu | Activation of α7 nicotinic acetylcholine receptors prevents monosodium iodoacetate-induced osteoarthritis in rats[END_REF][START_REF] Shaw | Janus kinase 2, an early target of alpha 7 nicotinic acetylcholine receptor-mediated neuroprotection against Abeta-(1-42) amyloid[END_REF]. Finally, we did not use Ach in vitro because of its short half-life and of its lack of specificity since Ach activates both nicotinic and muscarinic also the muscarinic receptors that might be expressed by chondrocytes.

We confirmed the anti-inflammatory and anti-catabolic effect of α7nAChR on human OA chondrocytes using a widely used and published specific α7nAChR agonist, PNU-282987 [START_REF] Grønlien | Distinct profiles of alpha7 nAChR positive allosteric modulation revealed by structurally diverse chemotypes[END_REF][START_REF] Hajós | The selective alpha7 nicotinic acetylcholine receptor agonist PNU-282987 [N-[(3R)-1-Azabicyclo[2.2.2]oct-3yl]-4-chlorobenzamide hydrochloride] enhances GABAergic synaptic activity in brain slices and restores auditory gating deficits in anesthetized rats[END_REF].

This result corroborates a previous study that showed the promotion of collagen synthesis by nicotine in human normal chondrocytes [START_REF] Gullahorn | Smoking and osteoarthritis: differential effect of nicotine on human chondrocyte glycosaminoglycan and collagen synthesis[END_REF]. However, the inflammatory and catabolic responses of human chondrocytes may depend on the expression level of the duplicate CHRFAM7A. We found that the expression level of the human-specific duplicate CHRFAM7A, which is responsible for a truncated protein of the α7 subunit, is associated with an enhanced catabolic profile (i.e., increased MMP expression). This duplicate contains exons 5-10 of the native gene of the alpha 7 subunit Chrna7, which transcribed most of the protein structure. However, the ligand binding site of this truncated protein is missing, leading to a decreased response to specific ligands, such as ACh or nicotine [START_REF] Chan | CHRFAM7A alters binding to the neuronal alpha-7 nicotinic acetylcholine receptor[END_REF]. The local expression of CHRFAM7A is suspected to be associated with

Accepted Article

This article is protected by copyright. All rights reserved some disease and inflammatory states, as demonstrated by the colon expression of CHRFAM7A in inflammatory bowel diseases [START_REF] Baird | Up-regulation of the humanspecific CHRFAM7A gene in inflammatory bowel disease[END_REF]. Interestingly, meta-analysis of observational studies have reported that smoking could be associated with less incident OA [START_REF] Kong | Association between smoking and risk of knee osteoarthritis: a systematic review and meta-analysis[END_REF]. Our study could give a rationale to these data. However, nicotine is just one of the compounds present in tobacco smoke and so the role of the other ones cannot be ruled out. Altogether, we still consider that even if nicotine could have some beneficial effects on cartilage, the multiple risks related to smoking far outweighs the conditional beneficial effect of nicotine.

Finally, we confirmed the involvement of α7nAChR in OA in vivo. We observed that Chrna7 -/- mice displayed more severe cartilage lesions than WT mice after MNX induction. This result confirms that mice lacking Chrna7 are associated with increased chondrocyte activation and cartilage degradation compared with WT mice and that α7nAChR has a protective role in OA.

This is the first study evaluating OA in Chrna7 -/-mice. These mice do not have behavioural and phenotypical difference compared to WT [START_REF] Yin | Chrna7 deficient mice manifest no consistent neuropsychiatric and behavioral phenotypes[END_REF]. However, we did not evaluate pain assessment. [START_REF] Teng | Nicotine Attenuates Osteoarthritis Pain and Matrix Metalloproteinase-9 Expression via the α7 Nicotinic Acetylcholine Receptor[END_REF][START_REF] Liu | Activation of α7 nicotinic acetylcholine receptors prevents monosodium iodoacetate-induced osteoarthritis in rats[END_REF]. They suggested that this effect was due to activation α7nAchR since the protective effect of nicotine disappears after injection of methyllycaconitine (MLA), an α7nAChR antagonist; however, the specificity of this chemical compound is questionable, and the administration of the combination of nicotine and MLA is somewhat artificial and potentially toxic to rats. In contrast to the more severe OA lesions caused by nicotine and MLA treatment in the MIA model, the exacerbated OA lesions in Chrna7 -/-MNX mice could be due to a loss of effect of endogen ACh on α7nAChR in Chrna7 -/-chondrocytes. A weakness of this study is that since we used a model of full invalidation of Chrna7 gene and not a cartilage specific one, we cannot definitely conclude whether cartilage genetic deletion is the only factor involved in the exacerbated OA lesions in Chrna7 -/-mice. Of note, it is possible that deletion of Chrna7 in bone also contributes to the observed MNX phenotype.

Our study has several strengths. First, we confirmed α7nAChR expression using various robust methods (i.e, α-bungarotoxin). Antibody-based techniques for the assessment of α7nAChR protein expression was not possible because of the unavailability of accurate and specific antibodies raised against the nicotinic subunit (41). We also evaluated for the first time the local expression of the specific human duplicate. Second, we used a mechanical model, which is much more representative of human OA than the highly inflammatory MIA model, far from human

Accepted Article

This article is protected by copyright. All rights reserved disease. Third, we used genetically modified mice instead of applying a combination of chemical compounds to mimic α7nAChR inhibition.

In conclusion, human and murine chondrocytes belong to the non-neuronal cholinergic system.

The alpha7 nicotinic ACh receptor is a crucial receptor able to limit the inflammatory and catabolic response of chondrocytes, and lacking this receptor can aggravate OA cartilage damage.

In human chondrocytes, expression of the gene duplicate CHRFAM7A is associated with a more marked catabolic profile, while its role in OA remains to be investigated. Thus, the non-neuronal

ACh/α7nAChR axis is involved in the pathophysiology of OA and could be a therapeutic target.

Accepted Article

This article is protected by copyright. All rights reserved

Accepted Article

This article is protected by copyright. All rights reserved A: 4X and 10 X representative Safranin O/fast green-stained sections of the right knee 8 weeks after sham surgery and MNX in WT and Chrna7 -/-mice .B: Weight of WT (n=9) and Chrna7 -/- (n=7) mice 8 weeks after sham surgery or MNX (not significant between groups). C: OARSI score (0-6, 6 being the worst) D: Synovitis score (0-6, 6 being the worst) and E : osteophyte score (0-3, 3 being the worst) between WT sham (n=6), Chrna7 -/-sham (n=5), WT MNX (n=9) and

Figure 1 :

 1 Figure 1: Non-neuronal cholinergic system expression in murine and human chondrocyte primary cultures and the in vitro effect of nicotine. A: Representative RT-PCR analysis of the cholinergic system expressed by murine chondrocytes (left) and human OA chondrocytes (right). We used water RNase-free conditions as a negative control for each primer (H20). B: qPCR expression of IL6 (left) and MMP3 (right) by WT murine primary culture of chondrocytes pretreated with 1, 10 or 100 µM nicotine for 15 minutes and stimulated with 10 ng/mL IL1β for 24 h (n=5 for WT, *p<0.05 using a non-parametric test). The results are normalized to the IL1β condition. C: ELISA results of IL6 (left) and MMP3 (right) in WT murine primary culture of chondrocytes pretreated with 1, 10 or 100 µM nicotine for 15 minutes and stimulated with 10 ng/mL IL1β for 24 h (n=5 for WT, *p<0.05 using a nonparametric test). The results are normalized to the IL1β condition. HPRT: hypoxanthine guanine phosphoribosyl transferase; ChAT: choline acetyltransferase; carAT: carnityl acetyltransferase; VAChT: vesicular acetylcholine transporter; AChE: acetylcholine esterase; BChE: butyrylcholine esterase; CHT1: choline transporter-1. IL6: interleukin-6; MMP3: metalloprotease-3, bp : base pair lenght.

Figure 2 :

 2 Figure 2: Nicotinic acetylcholine receptor expression by murine and human OA chondrocytes A: qPCR results of nicotinic acetylcholine receptor (nAChR) subunit expression by murine chondrocytes (left) and human OA chondrocytes (right). B: Representative and semi-quantitative α-bungarotoxin binding to α7nAChR on the surface of murine chondrocytes in primary culture.Cells were incubated with biotinylated α-bungarotoxin at 20 µg/mL overnight at 4°C and analyzed using a streptavidin and DAB protocol and a standard microscope, 20x. The control represents the same experiment without α-bungarotoxin incubation. Two independent experiments (one in red and one in green) were performed and each experiment was analysed in 4 different area of 91,000 µm 2 C: Representative and semi-quantitative α-bungarotoxin-binding α7nicotinic receptors on the surface of human OA chondrocytes in primary culture. Cells were incubated with biotinylated αbungarotoxin at 5 µg/mL. Two independent experiments (one in red and one in green) were performed and each experiment was analysed in 4 different area of 91,000 µm 2 .

Figure 3 :

 3 Figure 3: Expression of native Chrna7 and the specific duplicate CHRFAM7A by human OA chondrocytes A: CHRNA7 is the human gene encoding the a7nAChR; it is localized on chromosome 15q13,14 and contains 10 exons. It is assumed that during evolution, a partial duplication of exons 5-10 occurred approximately 1.6 MB upstream on the same chromosome. Thereafter, 4 completely novel exons (A-D) were rearranged in place of the missing exons 1-4 of human a7nAChR, and the 4 new exons of this FAM7 sequence created a new open reading frame to form a new hybrid gene called CHRFAM7A. Surprisingly, the "new" exons (A-D) were derived, in part, from an unrelated sequence, itself originating, in part, from human chromosome 3. Red represents the sequence amplified by the two primers. B: qPCR results of CHRNA7 (native) and CHRFAM7A expression by human OA chondrocytes (n=23, *p<0.05). C: Spearman correlation between CHRFAM7A duplicate mRNA expression and basal MMP3 and MMP13 mRNA expression (n=23, *p<0.05). D: MMP3 and MMP13 mRNA expression by human OA chondrocytes compared with CHRFAM7A mRNA expression by human OA chondrocytes (n=23, *p<0.05).

Figure 4 :

 4 Figure 4: The nicotine effect is mediated by α7nAchR A: qPCR expression and B: ELISA production of IL6 (left) and MMP3 (right) by Chrna7 -/-murine primary culture of chondrocytes pretreated with 1, 10 or 100 µM nicotine for 15 minutes and stimulated with 10 ng/mL IL1β for 24 h (n=4 for qPCR and n=5 for ELISA, *p<0.05 using a nonparametric test). The results are normalized to the IL1β condition. NS: Not significant. C: Calcium influx assay in primary WT and D: Chrna7 -/-chondrocyte culture. Calcium signaling was investigated by the use of a fluorescent Ca2+ indicator in murine WT (n=2) and Chrna7 -/-(n=2) chondrocytes treated with buffer alone, Ionomycin 10 µM and nicotine 10 µM. Analyzes were performed using a microscope Olympus IX83 inverted. E: Representative panel of the mean fluorescence measured immediately after buffer alone, Ionomycin 10 µM and nicotine 10 µM in WT and Chrna7 -/-chondrocytes. Nicotine did not elicit a transient increase in intracellular calcium. F: qPCR results of IL6, MMP3 and MMP13 expression by human OA chondrocytes stimulated with 10 ng/mL IL1β with or without a specific a7nAChR agonist, PNU 282987, at 100 µM (n=3).

Figure 5 :

 5 Figure 5: The lack of Chrna7 is associated with a more severe post-traumatic OA cartilage damage.

All rights reserved Human OA chondrocyte experiments. Human

  OA chondrocytes were isolated from the cartilage of patients undergoing joint arthroplasty for OA at Saint-Antoine Hospital or at the Clinique des

	Maussins-Ramsay (Paris, France), as previously described (19). Informed consent for the use of
	tissue was obtained from each patient before surgery. Experiments using human samples were
	approved by a French Institutional Review Board (Comité de Protection des Personnes, Paris Ile
	de France 5 and Commission Nationale de l'Informatique et des Libertés). Briefly, cartilage was
	dissected in small pieces, rinsed in saline buffer and digested with collagenases I and II (Roche,
	Penzberg, Germany). Chondrocytes were then seeded in plates at 20 × 10 3 cells per cm 2 and
	expanded in DMEM with 4.5 g/L glucose P/S/Glu containing 10% FCS. At confluence,
	chondrocytes were starved for 24 hours in DMEM with 4.5 g/L glucose P/S/Glu with 0.1% BSA
	and then stimulated for 24 h with IL-1β (10 ng/mL) (PeproTech, Rocky Hill, NJ, USA) with or
	without pretreatment with a specific 7 nicotinic receptor agonist, PNU-282987 (Abcam,
	Cambridge, U.K.), at 100 µM

  Liu et al. and Teng et al. previously reported that systemic administration of nicotine reduced monoiodoacetate (MIA) OA lesions

Acknowledgements: The authors acknowledge Professor Alain Sautet (MD) and Doctor

Francois-Paul Ehkirch (MD) for providing human knee joints, Pradeep Sacitharan (PhD) and Professor Pierre Aucouturier (MD, PhD) for their advice and helpful discussions, Kristell Wanderhick for her help in histological staining, Pierre Galichon and Dominique Prié for their advices on the calcium assay, and Tatiana Ledent for providing the animal facilities.

This work was supported by the French Society of Rheumatology (SFR) and the Assistance Publique -Hôpitaux de Paris (AP-HP). Laure Sudre was funded by the network Research of Osteoarthritis Diseases (ROAD).


Tahani Mtar
email: thani.mtar@gmail.com

Radhouane Fekih-Salem
email: radhouene.fekihsalem@isima.rnu.tn

Tewfik Sari
email: tewfik.sari@inrae.fr

INTERSPECIFIC DENSITY-DEPENDENT MODEL OF PREDATOR-PREY RELATIONSHIP IN THE CHEMOSTAT

Keywords: Bifurcation, chemostat, coexistence, predator-prey relationship, operating diagram. 10,

The objective of this study is to analyze a model of competition for one resource in the chemostat with general inter-specific density-dependent growth rates, taking into account the predator-prey relationship. This relationship is characterized by the fact that the prey species promotes the growth of the predator species which in turn inhibits the growth of the first species. The model is a three-dimensional system of ordinary differential equations. With the same dilution rates, the model can be reduced to a planar system where the two models have the same local and even global behavior. The existence and stability conditions of all steady states of the reduced model in the plane are determined according to the operating parameters. Using the nullcline method, we present a geometric characterization of the existence and stability of all equilibria showing the multiplicity of coexistence steady states. The bifurcation diagrams illustrate that the steady states can appear or disappear only through saddle-node or transcritical bifurcations. Moreover, the operating diagrams describe the asymptotic behavior of this system by varying the control parameters and show the effect of the inhibition of predation on the emergence of the bistability region and the reduction until the disappearance of the coexistence region by increasing this inhibition parameter.

Introduction

The culture of microorganisms occurs in a laboratory apparatus called "chemostat", which was invented in 1950. It is often used in various scientific fields such as ecology, biology, genetics, etc. It is also used as a model of the wastewater treatment process. We recall the classical chemostat model when two species x 1 and x 2 are in competition for a single substrate S:

       Ṡ = D(S in -S) - 1 γ 1 f 1 (S)x 1 - 1 γ 2 f 2 (S)x 2 , ẋ1 = (f 1 (S) -D)x 1 , ẋ2 = (f 2 (S) -D)x 2 , (1.1)
where S(t) is the concentration of the substrate at time t; x 1 (t) and x 2 (t) the concentrations of the species at time t; D and S in are, respectively, the dilution rate and the concentration of the substrate in the feed device; f 1 (S) and f 2 (S) are the growth functions and γ 1 , γ 2 the yield constants, which can be easily normalized to 1, using the change of variable x i /γ i → x i , i = 1, 2. Since at steady state, we cannot have, in general, f 1 (S) = f 2 (S) = D, the model (1.1) predicts that at most one species can survive the competition, the one with the lowest break-even concentration. This classical result is called the Competitive Exclusion Principle (CEP) and which has been corroborated by the experiences of Hansen and Hubbell [START_REF] Hansen | Single-nutrient microbial competition : Qualitative agreement between experimental and theoretically forecast outcomes[END_REF]. But, the CEP contradicts the biodiversity that is observed in microbial ecosystems as well as in waste-water treatment processes and bioreactors [START_REF] Scheffer | Why plankton communities have no equilibrium: solutions to the paradox[END_REF][START_REF] Schmidt | Characterization of a three bacte-22 REFERENCES ria mixed culture in a chemostat: evaluation and application of a quantitative Terminal-Restriction Fragment Length Polymorphism (T-RFLP) analysis for absolute and species specific cell enumeration[END_REF]. Thus, the classical chemostat model is unable to reproduce reality even qualitatively and new hypotheses should be considered in order to reconcile the theory and the experimental results. Several mechanisms of coexistence in the chemostat were considered in the literature, such as the intra-and interspecific competition between the species [START_REF] Abdellatif | Competition for a single resource and coexistence of several species in the chemostat[END_REF][START_REF] Borsali | Contribution to the study of the effect of the interspecificity on a two nutrients competition model[END_REF][START_REF] Leenheer | Crowding effects promote coexistence in the chemostat[END_REF][START_REF] Fekih-Salem | Sur un modle de comptition et de coexistence dans le chmostat[END_REF][START_REF] Wolkowicz | Direct interference on competition in a chemostat[END_REF], the flocculation of the species [START_REF] Fekih-Salem | Extensions of the chemostat model with flocculation[END_REF][START_REF] Fekih-Salem | Emergence of coexistence and limit cycles in the chemostat model with flocculation for a general class of functional responses[END_REF][START_REF] Fekih-Salem | Properties of the chemostat model with aggregated biomass and distinct removal rates[END_REF][START_REF] Haegeman | How flocculation can explain coexistence in the chemostat[END_REF], and the densitydependence of the growth functions [6-8, 10, 21]. Mathematical study of several chemostat-like models can be found in the monograph [START_REF] Smith | The Theory of the Chemostat: Dynamics of Microbial Competition[END_REF], see also [START_REF] Harmand | The Chemostat: Mathematical Theory of Microorganism Cultures[END_REF].

On the other hand, some studies have shown that species can coexist considering predator-prey interactions in the chemostat, [START_REF] Boer | Food chain dynamics in the chemostat[END_REF][START_REF] Li | Simple food chain in a chemostat with distinct removal rates[END_REF][START_REF] Nasrin | Three species food web in a chemostat[END_REF]. For instance, in [START_REF] Boer | Food chain dynamics in the chemostat[END_REF][START_REF] Li | Simple food chain in a chemostat with distinct removal rates[END_REF], a predator-prey model in the chemostat has been analyzed with Holling's type-II functional responses and distinct removal rates, where the model presents a rich set of behaviors with coexistence, transcritical and Hopf bifurcations with occurrence of stable limit cycles. Li and Kuang [START_REF] Li | Simple food chain in a chemostat with distinct removal rates[END_REF] show that the behavior of the system can be affected according to the death rate of species where the system can present oscillatory solutions. In [START_REF] Nasrin | Three species food web in a chemostat[END_REF], authors have considered the case of more complex food webs where two preys consume a single nutrient and the predator consumes only the preys. With the general monotone response functions and the same dilution rates, it was shown that the coexistence steady state is stable as long as it exists.

In this context, we will consider an interspecific density-dependent model of two species competing for a single nutrient in the chemostat where the removal rates of the substrate and the biomass are equal. This model introduced by El Hajji et al. Interspecific density-dependent model 3 [START_REF] Hajji | Association between competition and obligate mutualism in a chemostat[END_REF] is written as follows:

   Ṡ = D(S in -S) -f 1 (S, x 2 )x 1 -f 2 (S, x 1 )x 2 , ẋ1 = (f 1 (S, x 2 ) -D)x 1 , ẋ2 = (f 2 (S, x 1 ) -D)x 2 , (1.2)
where x 1 (t) and x 2 (t) are, respectively, the concentrations of prey and predator at time t. We assume that f 1 depends only on S and x 2 , and f 2 depends only on S and x 1 , and are increasing in the variable S. To demonstrate how interspecific interference can explain the coexistence of competing species in chemostat, El Hajji [START_REF] Hajji | How can inter-specific interferences explain coexistence or confirm the competitive exclusion principle in a chemostat?[END_REF], considered three types of interspecific interactions: inhibition, mutualism and food web (or predator-prey relationship). The inhibition is characterized by the fact that each species inhibits the growth of the other species. The mutualism is characterized by the fact that each species promotes the growth of the other species. The predator-prey relationship is characterized by the fact that one species (the prey) promotes the growth of the other species (the predator) which in turn inhibits the growth of the prey species.

El Hajji and Rapaport [START_REF] Hajji | How mutual inhibition confirms competitive exclusion principle[END_REF] have considered the first type where f 1 is decreasing in x 2 and f 2 is decreasing in x 1 . In [START_REF] Hajji | How can inter-specific interferences explain coexistence or confirm the competitive exclusion principle in a chemostat?[END_REF][START_REF] Hajji | How mutual inhibition confirms competitive exclusion principle[END_REF], it was proved that the coexistence of two species is impossible which confirms the CEP. Fekih-Salem et al. [START_REF] Fekih-Salem | A density-dependent model of competition for one resource in the chemostat[END_REF] have provided an extension of the results in [START_REF] Hajji | How mutual inhibition confirms competitive exclusion principle[END_REF] to general intra-and interspecific density-dependent growth rates f i which depends on both species x 1 and x 2 , not only on x j where j = i, and distinct removal rates.

The second type is the obligate mutualistic relationship where the function f i is increasing in the variable x j . Using general monotonic functional responses, the model exhibits the coexistence where the positive steady state is not necessarily unique [START_REF] Hajji | Association between competition and obligate mutualism in a chemostat[END_REF].

Here, we will give a particular attention to the third type of interaction which describes the predator-prey relationship where the function f 1 is decreasing in the variable x 2 while the function f 2 is increasing in the variable x 1 . In this case, El Hajji [START_REF] Hajji | How can inter-specific interferences explain coexistence or confirm the competitive exclusion principle in a chemostat?[END_REF] has shown that system (1.2) may exhibit the coexistence and the bistability with a multiplicity of positive steady states.

Our main objective in this work is to give a quite comprehensive analysis of model (1.2) in the case of the predator-prey relationship. The existence and stability conditions of all steady states are provided according to the control parameters S in and D. Using the nullcline method, we present a geometric characterization that describes all steady states of the model and determines the nature of bifurcation and the emergence of coexistence steady states. Varying the operating parameter S in , the bifurcation diagrams show all steady states that can appear or disappear only through saddle-node or transcritical bifurcations. All these results have not been done in [START_REF] Hajji | How can inter-specific interferences explain coexistence or confirm the competitive exclusion principle in a chemostat?[END_REF]. Last but not least, we analyze the operating diagram to determine the behavior of system with respect to the parameters S in and D, and we study the effect of the inhibition of predation on the behavior of system in this operating plane (S in , D). These questions were not considered in [START_REF] Hajji | How can inter-specific interferences explain coexistence or confirm the competitive exclusion principle in a chemostat?[END_REF]. This paper is organized as follows. In the next section, we first introduce the general assumptions on the growth functions and we provide a preliminary result on positivity and boundedness of solutions. Then, we define the reduced model in the plane from the three-dimensional system (1.2). In Section 3, we analyze the reduced model by determining the conditions of existence and stability of steady states Then, we give the result of global stability of the three-dimensional system (1.2). In Section 4, the operating diagrams describe the existence and the stability of each steady state according to control parameters D and S in . Considering specific growth rates, numerical simulations are presented in Section 5. The conclusion is given in the last Section 6. All proofs of propositions and lemmas are reported in Appendix A. With specific growth rates satisfying the general assumptions, the number of positive extrema of some functions is given in Appendix B. All the parameter values used in simulations are provided in Appendix C.

Assumptions and reduced model

In what follows, we consider model (1.2). For i = 1, 2, j = 1, 2, i = j, we first make the following general assumptions on the growth rates f i .

(H0) f i : R 2 + -→ R + is continuously differentiable. (H1) f i (0, x j ) = 0, for all x j ≥ 0.

(H2)

∂f i ∂S (S, x j ) > 0, for all S ≥ 0, x 1 > 0 and x 2 ≥ 0.

(H3) ∂f 1 ∂x 2 (S, x 2 ) < 0 and ∂f 2 ∂x 1 (S, x 1 ) > 0, for all S > 0, x 1 ≥ 0 and x 2 ≥ 0.

(H4) f 2 (S, 0) = 0, for all S > 0.

Assumption (H1) means that the substrate is necessary for the growth of the two species. Assumption (H2) means that the growth rate of each species increases with the concentration of the substrate. Assumption (H3) means that the growth of the first species x 1 is inhibited by the second species x 2 , while the growth of second species x 2 increases with the presence of first species x 1 . Assumption (H4) means that the prey species x 1 is necessary for the growth of the predator species x 2 . These assumptions are satisfied by the following growth rates:

f 1 (S, x 2 ) = m 1 S K 1 + S 1 1 + x 2 /L 1 , f 2 (S, x 1 ) = m 2 S K 2 + S x 1 L 2 + x 1 , (2.1)
where m 1 , m 2 denote the maximum growth rates; K 1 , K 2 and L 2 denote the Michaelis-Menten constants; L 1 represents the inhibition factor due to x 2 for the growth of the species x 1 . In order to show that our model (1.2) preserves the biological meaning, we will show the following property.

Proposition 2.1. Assume that Hypotheses (H0), (H1) and (H4) hold. For any nonnegative initial condition, the solution of system (1.2) exists for all t ≥ 0, re-10, 2020 0:15 WSPC/INSTRUCTION FILE MFS˙IJB˙0909˙WC

Interspecific density-dependent model 5 mains nonnegative and is bounded. In addition, the set

Ω = (S, x 1 , x 2 ) ∈ R 3 + : S + x 1 + x 2 = S in ,
is positively invariant and is a global attractor for the dynamics (1.2).

In what follows, it is convenient to use the abbreviation LES for Locally Exponentially Stable steady state and GAS for Globally Asymptotically Stable steady state, in the positive cone. In the following, we prove that the three-dimensional system (1.2) can be reduced to a two-dimensional system. Let z = S + x 1 + x 2 be the total mass density. Writing in the variables (z, x 1 , x 2 ), system (1.2) becomes

   ż = D(S in -z), ẋ1 = (f 1 (z -x 1 -x 2 , x 2 ) -D)x 1 , ẋ2 = (f 2 (z -x 1 -x 2 , x 1 ) -D)x 2 .
(2.2) Let (z(t), x 1 (t), x 2 (t)) be a solution of (2.2). From the first equation of (2.2), we deduce that z(t) = S in + (z(0) -S in )e -Dt . Thus, (x 1 (t), x 2 (t)) is a solution of the following asymptotically autonomous system:

ẋ1 = f 1 S in + (z(0) -S in )e -Dt -x 1 -x 2 , x 2 -D x 1 , ẋ2 = f 2 S in + (z(0) -S in )e -Dt -x 1 -x 2 , x 1 -D x 2 .
(2.3)

Then, system (2.3) converges toward the following autonomous system:

ẋ1 = (f 1 (S in -x 1 -x 2 , x 2 ) -D)x 1 , ẋ2 = (f 2 (S in -x 1 -x 2 , x 1 ) -D)x 2 . (2.4)
This is called the reduced model since it is simply the restriction of (2.3) to the invariant plane defined by z(t) = S in , for all t ≥ 0. Using Thieme's result [START_REF] Thieme | Convergence results and a Poincar-Bendixson trichotomy for asymptotically autonomous differential equations[END_REF], one can prove that the asymptotic behavior of the three-dimensional system (1.2) can be deduced from the asymptotic behavior of the two-dimensional system (2.4).

In the following, we focus on the study of the reduced system (2.4). This system is defined on the set

M := (x 1 , x 2 ) ∈ R 2 + : x 1 + x 2 ≤ S in ,
which is positively invariant for system (2.4). Indeed, for i = 1, 2, if x i = 0 then ẋi = 0. As a result, the axis x i = 0 is invariant under system (2.4) because the function

t -→ (0, x 2 (t)), (resp. t -→ (x 1 (t), 0)),
is a solution of (2.4). By uniqueness of solutions, the axis x i = 0 cannot be reached in finite time by trajectories for x i (0) > 0. Furthermore, we have

(x 1 + x 2 )(t 1 ) = S in , for some t 1 ≥ 0 =⇒ ẋi (t 1 ) = -Dx i (t 1 ) ≤ 0.
3. Analysis of the reduced and the three-dimensional models

The steady states of (2.4) are given by the solutions of the following system:

0 = (f 1 (S in -x 1 -x 2 , x 2 ) -D)x 1 , 0 = (f 2 (S in -x 1 -x 2 , x 1 ) -D)x 2 . (3.1)
From (H4), we deduce that system (2.4) cannot have a steady state where the first population is extinct while the second species survives. Otherwise, if x 1 = 0 one has x 2 = 0, then the steady state E 0 = (0, 0), called the washout, always exists. Assume that Hypotheses (H0) to (H2) hold. For all 0 D < f 1 (+∞, 0), we define the usual break-even concentration by

S = λ(D), (3.2)
which is the unique solution of equation

f 1 (S, 0) = D. If D ≥ f 1 (+∞, 0),
S in > λ(D). (3.3)
In the following, we determine the existence condition of the positive steady state

E * = (x * 1 , x * 2 ) of system (2.4)
, where both populations survive. The components x 1 = x * 1 and x 2 = x * 2 of E * must be the solutions of (3.1) with x 1 > 0 and x 2 > 0. Therefore, x * 1 and x * 2 are the solutions of the set of the following equations

f 1 (S in -x 1 -x 2 , x 2 ) = D, f 2 (S in -x 1 -x 2 , x 1 ) = D. (3.4)
That is to say, (3.4) must have a positive solution (x * 1 , x * 2 ) in

M o := (x 1 , x 2 ) ∈ R * + 2 : x 1 + x 2 < S in ,
which is the interior of M. To solve (3.4) in M o , we state the following results: Lemma 3.2. Assume that hypotheses (H0) to (H3) and condition (3.3) hold. The equation

f 1 (S in -x 1 -x 2 , x 2 ) = D defines a smooth decreasing function F 1 : [0, x1 ] -→ [0, x2 ] x 1 -→ F 1 (x 1 ) = x 2 , such that F 1 (x 1 ) = 0, F 1 (0) = x2 and -1 < F 1 (x 1 ) < 0 for all x 1 ∈ [0, x1 ],
where

x 2 = F 1 (x 1 ) is the unique solution of the equation f 1 (S in -x 1 -x 2 , x 2 ) = D and x2 is the unique solution of the equation f 1 (S in -x 2 , x 2 ) = D. Furthermore, the graph γ 1 of F 1 lies in M o (see Fig. 1 (b)). More precisely, (x 1 , F 1 (x 1 )) ∈ M o for all x 1 ∈]0, x1 [.
To determine the function defined by the equation f 2 (S in -x 1 -x 2 , x 1 ) = D, we will have need to define the solutions of the equation f 2 (S in -x 1 , x 1 ) = D. MFS˙IJB˙0909˙WC

Interspecific density-dependent model 7 Lemma 3.3. Assume that assumptions (H0) to (H4) hold. The equation

f 2 (S in -x 1 , x 1 ) = D (3.5)
has a solution in [0, S in ] if and only if,

max x1∈[0,Sin] f 2 (S in -x 1 , x 1 ) ≥ D. (3.6)
Generically, we have an even number of solutions in [0, S in ]. For simplicity, we add the following assumption which is satisfied by the specific growth rates (2.1).

(H5) Equation (3.5) has at most two solutions x 1 1 and x 2 1 in [0, S in ] (see Fig. 1(a)). With the specific growth rates (2.1), this property is satisfied (see Appendix B). The general case where the function x 1 → f 2 (S in -x 1 , x 1 ) is multimodal can be treated similarly, without added difficulty. In this particular case, we obtain the following result: Lemma 3.4. Assume that hypotheses (H0) to (H5) and condition (3.6) hold. The equation f 2 (S in -x 1 -x 2 , x 1 ) = D defines a function

F 2 : x 1 1 , x 2 1 -→ [0, S in [ x 1 -→ F 2 (x 1 ) = x 2 , such that F 2 x 1 1 = F 2 x 2 1 = 0 and -1 < F 2 (x 1 ) for all x 1 ∈ x 1 1 , x 2 1
, where andx 1 1 and x 2 1 are the solutions of (3.5). Furthermore, the graph

x 2 = F 2 (x 1 ) is the unique solution of the equation f 2 (S in -x 1 -x 2 , x 1 ) = D,
γ 2 of F 2 lies in M o (see Fig. 1(c)). More precisely, (x 1 , F 2 (x 1 )) ∈ M o for all x 1 ∈ x 1 1 , x 2 1 .
Note that the function F 2 has a unique extremum in x 1 1 , x 2 1 with the specific growth functions (2.1) (see Fig. 1(c) and Appendix B). One can see that x1 , x 1 1 and x 2 1 are simply the solutions of the following equations

f 1 (S in -x1 , 0) = D, f 2 S in -x 1 1 , x 1 1 = D, f 2 S in -x 2 1 , x 2 1 = D.
These quantities represent the coordinates of the intersections of the curves γ 1 and γ 2 with the x 2 = 0 axis. Their relative positions play a major role in the behavior of system (2.4). In fact, there exists three cases that must be distinguished:

Case 1 : x1 < x 1 1 < x 2 1 , Case 2 : x 1 1 < x1 < x 2 1 , Case 3 : x 1 1 < x 2 1 < x1 . (3.7)
When hypotheses (H0) to (H5) and condition (3.6) hold, we have the following equivalences (see Fig. 1(a)):

f 2 (λ(D), S in -λ(D)) < D ⇔ Case 1 or Case 3 f 2 (λ(D), S in -λ(D)) > D ⇔ Case 2. (3.8)
As a consequence of the previous lemmas, we obtain the following result which gives the condition of existence of positive steady states. ) is a positive solution of equations

x 2 = F 1 (x 1 ) and x 2 = F 2 (x 1 ). ( 3.9)
The following proposition provides the number of positive steady states of (2.4).

Proposition 3.6. Assume that hypotheses (H0) to (H5), and conditions (3.3) and (3.6) hold.

(1) When Case 1 holds, then there is no positive steady state.

(2) When Case 2 holds, then there exists at least one positive steady state. Generically, the system has an odd number of positive steady states. (3) When Case 3 holds, then generically system (2.4) has no positive steady state or an even number of positive steady states.

(a) The following result provides the local stability condition of each steady state of system (2.4).

x1 x2 γ1 γ2 δ x1 x 1 1 x 2 1 E1 • E0 • (b) γ2 γ1 x1 x2 δ x1 x 1 1 x 2 1 E1 • E0 • E * 1 • (c) γ2 γ1 x 1 1 x 2 1 E * 1 • E * 2 • E1 • E0 • x1 x2 δ

10, 2020 0:15 WSPC/INSTRUCTION FILE MFS˙IJB˙0909˙WC

Interspecific density-dependent model 9 Proposition 3.7. Assume that hypotheses (H0) to (H5) hold, we have

(1) E 0 = (0, 0) is LES if and only if S in < λ(D). (2) E 1 = (x 1 , 0) is LES if and only if f 2 (λ(D), S in -λ(D)) < D. (3) E * = (x * 1 , x * 2 ) is LES if and only if F 1 (x * 1 ) < F 2 (x * 1
). The next proposition determines the conditions of global asymptotic stability of each steady state of (2.4). In what follows, we will determine the global asymptotic stability of each steady state of three-dimensional system (1.2), which is deduced from that of reduced model (2.4). Notice that to the steady states of (1.2)

E 0 = (S in , 0, 0), E 1 = (S in -x1 , x1 , 0) , E * = (S in -x * 1 -x * 2 , x * 1 , x *
2 ) correspond the steady states of (2.4), respectively,

E 0 = (0, 0), E 1 = (x 1 , 0) , E * = (x * 1 , x * 2 )
. Using the convergence theorem [START_REF] Smith | The Theory of the Chemostat: Dynamics of Microbial Competition[END_REF] (or Thieme's result [START_REF] Thieme | Convergence results and a Poincar-Bendixson trichotomy for asymptotically autonomous differential equations[END_REF]), we deduce that the asymptotic behavior of the solution of three order system (1.2) is the same as the asymptotic behavior of reduced system (2.4). Now, the main result can be stated. Theorem 3.9. Assume that hypotheses (H0) to (H5) hold.

(1) E * is GAS in the interior of Ω if and only if f 2 (λ(D), S in -λ(D)) > D and (3.9) has a unique solution. The necessary and sufficient conditions of existence, local and global stability of all steady states of (1.2) are summarized in Table 1.

Table 1. Necessary and sufficient conditions of existence and stability of steady states of system (1.2).

Existence

Local stability Global stability

E 0 always exists S in < λ(D) S in ≤ λ(D) E 1 S in > λ(D) f 2 (λ(D), S in -λ(D)) < D f 2 (λ(D), S in -λ(D)) ≤ D, (3.
9) has no solution E * (3.9) has a solution

F 2 x * 1 > F 1 x * 1 f 2 (λ(D), S in -λ(D)) > D, (3.
9) has a unique solution

Operating diagram

The operating diagram shows how the system behaves when we vary the two control parameters S in and D. It will be used to describe the behavior of system (1.2) in different regions of the operating plan. All other biological parameters are those of the growth functions and are fixed since they depend on the nature of the organisms and the substrate introduced into the bioreactor. All the biological parameter values used are provided in Table 3 by considering the specific growth functions (2.1). Note that the operating diagrams of the interspecific density-dependent model (1.2) have not been studied in the existing literature. Let Υ 1 be the curve of equation

D = f 1 (S in , 0), Υ 1 := {(S in , D) : D = f 1 (S in , 0)} .
Recall that x1 = S in -λ(D) and x i 1 , i = 1, 2 are the solutions of equation (3.5) defined by f 2 (S in -x 1 , x 1 ) = D. Therefore, x1 and x i 1 depend on the operating parameters S in and D. We denote them by x1 (S in , D) and

x i 1 (S in , D), i = 1, 2. Define the curve Υ 2 of equation x1 (S in , D) = x 2 1 (S in , D) by the set Υ 2 := (S in , D) : x1 (S in , D) = x 2
1 (S in , D) . Let Υ 3 be the curve in the plan (S in , D) which corresponds to a saddle-node bifurcation with the appearance of two positive steady states. In this limit case, the curve γ 1 of the function x 1 → F 1 (x 1 ) is tangent to the curve γ 2 of the function x 1 → F 2 (x 1 ). Consequently, the equation F 1 (x 1 ) = F 2 (x 1 ) can have a solution x 1 on [0, S in ] which depends on S in and D, where this solution x 1 (S in , D) must satisfy a second equation F 1 (x 1 ) = F 2 (x 1 ) whose solution can define the curve Υ 3 . Thus, we define the set which is a curve in the generic case (see Fig. 4). Table 2 shows the existence and local stability of steady states E 0 , E 1 , E *

Υ 3 := {(S in , D) : F 1 (x 1 ) = F 2 (x 1 ) and F 1 (x 1 ) = F 2 (x 1 )} ,
1 and E * 2 in four regions labeled as J k , k = 0, . . . , 3, of the operating diagram shown in Fig. 4. In Table 2, the letter S (resp. U) means stable (resp. unstable) steady state. Absence of letter means that the corresponding steady state does not exist.

Table 2. Existence and local stability of steady states according to regions in the operating diagram of Fig. 4. Fig. 4 illustrates the operating diagram of (1.2) where the curves Υ i , i = 1, 2, 3 separate the plane (S in , D) in at most four regions J k , k = 0, . . . , 3. All regions appear in Figs. 4 (b) and 5 (b). The region J 3 (resp. J 2 ) is empty in case (c) in Fig.

Region E 0 E 1 E * 1 E * 2 (S in , D) ∈ J 0 S (S in , D) ∈ J 1 U S (S in , D) ∈ J 2 U S S U (S in , D) ∈ J 3 U U S (

(resp. (a)).

The transition from the region J 0 to the region J 1 by the curve Υ 1 corresponds to a transcritical bifurcation making the steady state E 0 unstable with the appearance of an LES steady state E 1 . The transition from the region J 1 to the region J 3 by the curve Υ 2 corresponds to a transcritical bifurcation making E 1 unstable with the appearance of the positive steady state E * which is LES. The transition from the region J 1 to the region J 2 by the curve Υ 3 corresponds to a saddle-node bifurcation with the appearance of two positive steady states E * 1 which is LES and E * 2 which is unstable. The transition from the region J 2 to the region J 3 by the curve Υ 2 corresponds to a transcritical bifurcation when the unstable steady state E * 2 disappears and E 1 becomes unstable. Finally, the operating diagram in Fig. 5 shows the effect of the parameter L 2 in the growth function (2.1) on the emergence of the bistability region J 2 (in yellow). When the parameter L 2 increases, the operating diagram illustrates how the coexistence region J 3 (in red) is reduced and even disappears. To maintain the coexistence of species, the ideal control parameter values of D and S in should be chosen in J 3 .

To expose all bifurcations already mentioned and the effect on the concentrations of species at steady state as the control parameter S in varies while D = 0.5, we study the one-parameter bifurcation diagram for system (1.2). Fig. 6 shows that the washout steady state E 0 loses its stability by a transcritical bifurcation with E 1 at the green diamond point for S in 0.44. Then, a saddle-node bifurcation occurs at cyan circle point for S in 3.29 between E * 1 and E * 2 . Finally, the steady state E 1 loses its stability by a transcritical bifurcation with E * 2 which disappears at the green diamond point for S in 9.35. • E1 E * 2 Fig. 6. One-parameter bifurcation diagrams when D = 0.5 showing the effect on the components x 1 and x 2 of all steady states as S in varies. Blue dashed curves correspond to unstable steady states and red solid curves to LES steady states. The green solid diamonds represent the transcritical bifurcations while the cyan solid circles represent the saddle-node bifurcations.

Numerical simulations

To validate the results of the previous section, we consider the growth functions given by (2.1) and the set of parameter values provided in Table 3 for plotting the operating diagram of Fig. 4. In this section, we plot the solutions of the reduced system (2.4) in the variables x 1 and x 2 . The phase portraits for system (2.4) illustrate the solutions in the x 1 -x 2 plane for several initial conditions. (a) 2 and E 0 are unstable when (S in , D) = (2.9, 0.3) is chosen in the bistability region J 2 . In this case, the numerical simulations show the bistability with two basins of attraction, which are separated by the stable manifold of saddle point E * 2 . One basin attracts the solutions to the positive steady state E * 1 and the other to the steady state E 1 . Fig. 7 (c) illustrates the case when (S in , D) = (10, 0.45) ∈ J 3 where there exists a unique positive steady state E * (0.381, 4.165), which is LES. Numerical simulations show the global convergence to E * for any positive initial condition.

x1 x2 E1 • E0 • (b) x1 x2 E1 • E0 • E * 1 • E * 2 • (c) x1 x2 E1 • E0 • E * •

Conclusion

In this work, we have analyzed mathematically and through numerical simulations an interspecific density-dependent model of two microbial species in a chemostat competing for a nutrient proposed by El Hajji [START_REF] Hajji | How can inter-specific interferences explain coexistence or confirm the competitive exclusion principle in a chemostat?[END_REF], taking into account the predatorprey relationship. With the same dilution rates, three order model (1.2) can be reduced to second order model (2.4). Thus, we have studied the reduced model (2.4) since the asymptotic behavior of each steady state of (1.2) is derived from the reduced model (2.4). Under general monotonic density-dependent growth functions, we give, in a first step, a complete analysis for the existence and stability of all steady states of model (2.4) by using the nullcline method. We have shown the multiplicity of positive steady states, where the number of positive steady states depends on the relative positions of the numbers x1 , x 1 1 and x 2 1 , namely, three cases of (3.7) must be distinguished. Furthermore, the stability of a positive steady state E * depends on the relative positions of the curves of the nullclines defined by (3.4).

On the other hand, we determined the operating diagram which depicts the asymptotic behavior of the system with respect to control parameters. It shows the occurrence of the bistability region with multiplicity of coexistence steady states that can bifurcate only through saddle-node bifurcations or transcritical bifurcations. To ensure the coexistence between the prey and predator species in the chemostat, the ideal control parameter values of D and S in should be chosen in the coexistence region.

The increase in the value of L 2 , that is, the decrease in the dynamics of predation reflects on the operating diagram the direct impact of the inhibition of predation on the prey. Indeed, a bistability region emerges while the coexistence is reduced and then disappears by increasing this inhibition parameter. When only S in varies, the one-parameter bifurcation diagram illustrates all bifurcations and shows the effect on the concentrations of prey and predator species at steady state. The numerical simulations show that the system may exhibit bistability with convergence either to the coexistence of prey and predator species or to the exclusion of the predator species when the growth functions are given by (2.1). function lies in M o (see Fig. 1(b)). By the implicit function theorem, the function F 1 is smooth and decreasing. Indeed, using (H2)-(H3), it follows that

-1 < F 1 (x 1 ) = - ∂f1 ∂x1 (x 1 , F 1 (x 1 )) ∂f1 ∂x2 (x 1 , F 1 (x 1 )) = ∂f1 ∂S -∂f1 ∂S + ∂f1 ∂x2 < 0.
Furthermore, if x 2 = 0 then

F 1 (x 1 ) = 0 ⇔ f 1 (S in -x 1 , 0) = D ⇔ x 1 = x1 = S in -λ(D),
and if x 1 = 0, then

F 1 (0) = x 2 ⇔ f 1 (S in -x 2 , x 2 ) = D.
From (H0) to (H3), we see that the function

x 2 -→ f 1 (S in -x 2 , x 2 )
is decreasing from f 1 (S in , 0) for x 2 = 0 to 0 for x 2 = S in . Consequently, the equation

f 1 (S in -x 2 , x 2 ) = D has a solution x2 ∈]0, S in [ if and only if f 1 (S in , 0) > D,
f 2 (S in -x 1 , x 1 ) ≥ D,
which is the condition (3.6).

Proof of Lemma 3.4. Assume that hypotheses (H0) to (H5) and condition (3.6) hold. From Lemma 3.3, the equation f 2 (S in -x 1 , x 1 ) = D has two solutions x 1 1 and

x 2 1 in [0, S in ]. Let us prove that (x 1 , F 2 (x 1 )) ∈ M o for all x 1 ∈ x 1 1 , x 2 1 .
Let l 1 be a fixed line defined by x 1 = b 1 that intersects the line δ defined by x 1 + x 2 = S in at point (b 1 , c 2 ) where c 2 = S in -b 1 . From hypotheses (H0) to (H2), we see that the function

x 2 -→ f 2 (S in -b 1 -x 2 , b 1 ) is decreasing from f 2 (S in -b 1 , b 1 ) for x 2 = 0 to 0 for x 2 = c 2 . Therefore, there exists a unique b 2 ∈ [0, c 2 [ such that f 2 (S in -b 1 -b 2 , b 1 ) = D if and only if f 2 (S in -b 1 , b 1 ) ≥ D, which is equivalent to b 1 ∈ x 1 1 , x 2
1 , using (H0) to (H5) and condition (3.6) (see Figs. 1(a,c)). Therefore, we have shown that for all b

1 ∈ x 1 1 , x 2 1 there exists a unique b 2 ∈ [0, c 2 ] such that f 2 (S in -b 1 -b 2 , b 1 ) = D.
-1 < F 2 (x 1 ) = - ∂f2 ∂x1 (x 1 , F 2 (x 1 )) ∂f2 ∂x2 (x 1 , F 2 (x 1 )) = -1 + ∂f2 ∂x1 ∂f2 ∂S .
In addition, when x 2 = 0, we have

F 2 (x 1 ) = 0 is equivalent to f 2 (S in -x 1 , x 1 ) = D, that is, x 1 = x 1 1 or x 1 = x 2 1
. This completes the proof. Proof of Proposition 3.7. At E 0 , the Jacobian matrix of reduced model (2.4) is written as follows:

J E0 = f 1 (S in , 0) -D 0 0 -D .
The eigenvalues are -D and f 1 (S in , 0) -D, which are negative if and only if S in < λ(D).

The Jacobian matrix at E 1 is given by

J E1 = -x 1 ∂f1 ∂S x1 (-∂f1 ∂S + ∂f1 ∂x2 ) 0 f 2 (S in -x1 , x1 ) -D .
The eigenvalues are on the diagonal. They are negative if and only if

f 2 (λ(D), S in -λ(D)) < D,
which is equivalent to Case 1 or Case 3 holds, thanks to (3.8) (see Fig. 1(a)).

The Jacobian matrix at a positive steady state E * is given by

J E * =   -x * 1 ∂f1 ∂S -x * 1 ∂f1 ∂S + ∂f1 ∂x2 -x * 2 ∂f2 ∂S + ∂f2 ∂x1 -x * 2 ∂f2 ∂S   . Since tr J E * = -x * 1 ∂f 1 ∂S + x * 2 ∂f 2 ∂S < 0,
then the real part of the eigenvalues of J E * are negative if and only if

det J E * = x * 1 x * 2 ∂f 2 ∂S ∂f 1 ∂x 2 + ∂f 2 ∂x 1 ∂f 1 ∂S - ∂f 2 ∂x 1 ∂f 1 ∂x 2 > 0.
Using the derivatives expressions of F 1 and F 2 , we can write

F 1 -F 2 = ∂f2 ∂S ∂f1 ∂x2 + ∂f2 ∂x1 ∂f1 ∂S -∂f2 ∂x1 ∂f1 ∂x2 ∂f2 ∂S ∂f1 ∂x2 -∂f1 ∂S .
As a consequence, one has

det J E * = (F 2 -F 1 )x * 1 x * 2 ∂f 2 ∂S ∂f 1 ∂S - ∂f 1 ∂x 2 .
From (H2) and (H3), it follows that

det J E * > 0 ⇔ F 2 (x * 1 ) > F 1 (x * 1 ) .
Proof of Proposition 3.8. In order to show that the solutions of (2.4) converge to one of the steady states of the system for all t ≥ 0, we apply the method used in [START_REF] Fekih-Salem | A density-dependent model of competition for one resource in the chemostat[END_REF][START_REF] Hirsch | Differential Equations, Dynamical Systems, and Linear Algebra[END_REF]. The curves γ 1 and γ 2 separate the set M o in at most four regions as shown in Fig. 3(a). We define the regions as follows:

I : ẋ1 < 0, ẋ2 < 0, II : ẋ1 < 0, ẋ2 > 0, III : ẋ1 > 0, ẋ2 < 0, IV : ẋ1 > 0, ẋ2 > 0.
Consider a trajectory (x 1 (t), x 2 (t)) ∈ M o starting at a positive initial condition (x 1 (0), x 2 (0)). Let us first prove that if f 2 (λ(D), S in -λ(D)) > D, that is, Case 2 holds (see equivalence (3.8)) and (3.9) has a unique solution, then E * is GAS in M o . In this case, the curves γ 1 and γ 2 have only one intersection and separate the set M o in four regions I, II, III and IV , see Fig. 3(a). Assume first that (x 1 (0), x 2 (0)) ∈ I. If (x 1 (t), x 2 (t)) remains in I for all t ≥ 0 then x 1 (t) and x 2 (t) are positive decreasing functions, with limits If (x 1 (t), x 2 (t)) leaves the region I and enters in III. Then (x 1 (t), x 2 (t)) can be remained in III for all t ≥ 0, where x 1 (t) is increasing and x 2 (t) is decreasing. Thus, the limits (A.4) exist and (x 1∞ , x 2∞ ) is a steady state of (2.4) which belongs to III. Since E * and E 0 are the only steady states in III and since E 0 attracts only solutions with x 1 (0) = 0, we conclude that (A.5) holds. On the other hand, if (x 1 (t), x 2 (t)) leaves the region III and enters in IV , then also (x 1 (t), x 2 (t)) can be remained in IV for all t ≥ 0, where x 1 (t) and x 2 (t) are increasing. Thus the limits (A.4) exist and (x 1∞ , x 2∞ ) is a steady state of (2.4) which belongs to IV . Since E * et E 1 are the only steady states in IV and since E 1 attracts only solutions with x 2 (0) = 0, we conclude that (A.5) holds. The solution (x 1 (t), x 2 (t)) can leave the region IV and enters in II. If (x 1 (t), x 2 (t)) remains in the region II for all t ≥ 0 then x 1 (t) is a positive decreasing function and x 2 (t) is a positive increasing function. Therefore, the limits (A.4) exist and (x 1∞ , x 2∞ ) is a steady state of (2.4) which belongs to II. Since E * et E 1 are the only steady states in II and since E 1 attracts only solutions with x 2 (0) = 0, we conclude that (A.5) holds. Finally, if (x 1 (t), x 2 (t)) leaves the region II and enters in I, then as shown previously it necessarily tends to E * and hence (A.5) holds. Therefore, E * is GAS in M o . Inversely, if E * is GAS in M o , then only the Case 2 or Case 3 can be satisfied, using Proposition 3.6. The Case 3 is impossible since in this case we have E 1 is LES, using Proposition 3.7. Thus, Case 2 holds. If there exist at least two positive steady states, then E * cannot be GAS in M o which is a contradiction. Hence, there exists a unique positive steady state.

Then, the equation F 2 (x 1 ) = 0 has at most two solutions on x 1 1 , x 2 1 . In other hand, from Lemma 3.4, the function F 2 is positive and continuous on x 1 1 , x 2 1 and F 2 (x 1 1 ) = F 2 (x 2 1 ) = 0. Therefore, F 2 must have a unique positive extremum on x 1 1 , x 2 1 (see Fig. 1(c)).

Appendix C. Parameters used in numerical simulations

All parameter values used in numerical simulations are provided in Table 3.

2 : 3 . 1 .

 231 we put λ(D) = +∞. The following result determines the existence of the steady state corresponding to the extinction of the second population x Proposition Assume that assumptions (H0)-(H2) hold. System (2.4) has the steady state E 1 = (x 1 , 0) where x1 = S in -λ(D). It exists if and only if,

Fig. 1 .

 1 Fig. 1. When the specific growth functions are given by (2.1): (a) number of solutions of equation (3.5), (b) definition of F 1 , (c) definition of F 2 where δ is the line defined by x 1 + x 2 = S in .

Fig. 2 .

 2 Fig. 2. Number of positive steady states:(a) Case 1: no intersection, (b) Case 2: an odd number of intersections, (c) Case 3: an even number of intersections. In all figures, we have chosen the red color for LES steady state and the blue color for unstable steady state.

Proposition 3 . 8 .

 38 Assume that hypotheses (H0) to (H5) hold. (1) E * is GAS in M o if and only if f 2 (λ(D), S in -λ(D)) > D and (3.9) has a unique solution. (2) E 1 is GAS in M o if and only if f 2 (λ(D), S in -λ(D)) ≤ D and (3.9) has no solution. (3) E 0 is GAS in M o if and only if S in ≤ λ(D).

Fig. 3 .

 3 Fig. 3. Global stability of steady states of (2.4):(a) Item 1 of Proposition 3.8, where E * is GAS, (b) Item 2 of Proposition 3.8, where Case 1 holds and E 1 is GAS, (c) Item 3 of Proposition 3.8, where E 0 is GAS.

  Mtar, R. Fekih-Salem, T. Sari(2) E 1 is GAS in the interior of Ω if and only if f 2 (λ(D), S in -λ(D)) ≤ D and (3.9) has no solution. (3) E 0 is GAS in Ω if and only if S in ≤ λ(D).

Fig. 4 .

 4 Fig. 4. (a) The curves Υ i , i = 1, 2, 3 when Υ 2 and Υ 3 intersect in S * in , D * (37.95, 2.49). (b) The corresponding operating diagram of (1.2). (c) Bistability and multiplicity of positive steady states when (S in , D) = (3.8, 0.5) ∈ J 2 .

Fig. 5 .

 5 Fig. 5. Effect of variation of the parameter L 2 on the operating diagram: (a) L 2 = 0.1, (b) L 2 = 1.5, (c) L 2 = 10.

  0

Fig. 7 .

 7 Fig. 7. Phase portraits for (2.4) according to regions in the operating diagram of Fig. 4: (a) Global convergence to E 1 for (S in , D) ∈ J 1 . (b) Bistability of E * 1 and E 1 for (S in , D) ∈ J 2 . (c) Global convergence to E * for (S in , D) ∈ J 3 .

Fig. 7 ( 1 ( 1 .

 711 Fig. 7 (a) illustrates the global convergence to E 1 which is GAS in M o when the system has no positive steady state where (S in , D) = (2.5, 0.3) ∈ J 1 . Fig. 7 (b) illustrates the multiplicity of positive steady states and the bistability of E 1 (2.655, 0) and E * 1 (1.111, 1.223) which are LES while E *2 and E 0 are unstable when (S in , D) = (2.9, 0.3) is chosen in the bistability region J 2 . In this case, the numerical simulations show the bistability with two basins of attraction, which are separated by the stable manifold of saddle point E * 2 . One basin attracts the solutions to the positive steady state E * 1 and the other to the steady state E 1 . Fig.7(c) illustrates the case when (S in , D) = (10, 0.45) ∈ J 3 where there exists a unique positive steady state

  that is, condition(3.3) holds. If such an x2 exists, then it is unique. Proof of Lemma 3.3. Let φ 2 be the function defined byφ 2 (x 1 ) := f 2 (S in -x 1 , x 1 ). (A.3) Under hypothesis (H0)-(H4), φ 2 is positive and continuous on [0, S in ], nonmonotonic and φ 2 (0) = φ 2 (S in ) = 0. Therefore, the equation f 2 (S in -x 1 , x 1 ) = D has a solution in [0, S in ] if and only if max x1∈[0,Sin]

17 F 2

 172 by b 2 = F 2 (b 1 ). The graph γ 2 of this function lies in M o . Using the implicit function theorem and assumptions (H2)-(H3), we obtain

lim t→+∞ x 1

 1 (t) = x 1∞ and lim t→+∞ x 2 (t) = x 2∞ . (A.4)Therefore, (x 1∞ , x 2∞ ) is a steady state of (2.4) which belongs to the closure I of the set I. Since E * is the only steady state in I, henceE * = (x 1∞ , x 2∞ ). (A.5)

  Proposition 3.5. Assume that hypotheses (H0) to (H5), and conditions (3.3) and (3.6) hold. A positive steady state E * = (x * 1 , x * 2 ) of (2.4) exists if and only if the curves γ 1 and γ 2 have a positive intersection, that is, (x * 1 , x * 2

Table 3 .

 3 Parameter values used for (2.4) when f 1 and f 2 are given by (2.1).

	Figure	m 1	K 1	L 1	m 2	K 2	L 2	S in	D
	Figs. 1 (a,c)				2.95	1.8	1.5	5	0.5
	Fig. 1 (b)	2.75	2	1.2				5	0.5
	Fig. 2 (a)	2.5	2	1	6.7	1	2	2	1
	Fig. 2 (b)	1.2	1.8	2.2	3	1.8	1	3	0.5
	Fig. 2 (c)	2.3	1.8	2.2	2.2	1.8	1	5	0.5
	Figs. 4(a,b), 5(b) 2.75	2	1.2 2.95	1.8	1.5 variable variable
	Fig. 4 (c)	2.75	2	1.2 2.95	1.8	1.5	3.8	0.5
	Fig. 5 (a)	2.75	2	1.2 2.95	1.8	0.1	variable variable
	Fig. 5 (c)	2.75	2	1.2 2.95	1.8	10	variable variable
	Fig. 6	2.75	2	1.2 2.95	1.8	1.5	variable	0.5
	Fig. 7 (a)	2.75	2	1.2 2.95	1.8	1.5	2.5	0.3
	Fig. 7 (b)	2.75	2	1.2 2.95	1.8	1.5	2.9	0.3
	Fig. 7 (c)	2.75	2	1.2 2.95	1.8	1.5	10	0.45

T. Mtar, R. Fekih-Salem, T. Sari

T. Mtar, R. Fekih-Salem, T. Sari

T. Mtar, R. Fekih-Salem, T. Sari

Acknowledgments

The first author thanks the Tunisian Ministry of Higher Education, Scientific Research and Technology for financial support and INRAE for hosting her in Montpellier during the preparation of this work. This work was supported by the Euro-Mediterranean research network TREASURE (http://www.inra.fr/treasure). The authors thank the anonymous reviewers for their constructive comments which have greatly improved this work. 10, 2020 0:15 WSPC/INSTRUCTION FILE MFS˙IJB˙0909˙WC

Appendix A. Proofs

In this Appendix, we give the proofs of all propositions and lemmas.

Proof of Proposition 2.1. From (H0), the vector field defined by (1.2) is continuously differentiable. Hence, the existence and the uniqueness of solution to initial value problems hold. Thanks to (H1), if S = 0 then Ṡ = DS in > 0. Thus, no trajectory can leave the positive octant R 3 + by crossing the boundary face S = 0. In addition, if x 1 = 0 and x 2 = 0 then ẋ1 = 0 and ẋ2 = 0 and the set Interspecific density-dependent model 15 is a solution of (1.2). By uniqueness of solutions, the set Γ 0 cannot be reached in finite time by trajectories for which x 1 (0) > 0 or x 2 (0) > 0. Furthermore, for i = 1, 2, if x i = 0 then ẋi = 0 and the set

2) because the function

is a solution of (1.2), thanks to (H4). By uniqueness of solutions, the set Γ i cannot be reached in finite time by trajectories for x i (0) > 0. Therefore, the solutions remain non-negative. Let z = S + x 1 + x 2 be the total mass density. The sum of the three equations of (1.2) yields

We deduce that 0 ≤ z(t) < max(z(0), S in ), for all t ≥ 0.

Therefore, the solutions of (1.2) are positively bounded and are defined for all t ≥ 0. From (A.1), it can be deduced that the set Ω is positively invariant and is a global attractor for (1.2).

Proof of Proposition 3.1. If x 2 = 0, from the first equation of (3.1), x 1 is the solution of the equation

From Hypotheses (H0) to (H2) and definition (3.2) of λ, we have

Proof of Lemma 3.2. Let l 1 be a fixed line defined by x 1 = b 1 that intersects the line δ defined by

From hypotheses (H0) to (H3), we see that the function A similar argument to previous case shows that E 1 is GAS in M o in this case (see Fig. 3(b)). Inversely, if E 1 is GAS in M o , then the positive steady state E * does not exist and E 1 is locally asymptotically stable. From the Jacobian matrix at E 1 (see the proof of Proposition 3.7), one has

The last item of the proposition can be treated with the same manner (see Fig. 3(c)).

Appendix B. The particular case for growth functions (2.1)

The following result determines the maximal number of solutions of equation (3.5) in the particular case of growth functions (2.1).

Lemma Appendix B.1. Consider the specific growth functions (2.1) and assume that (3.6) holds. Then, equation (3.5) has at most two solutions x 1 1 and

Proof. When the growth functions are given by (2.1), from definition (A.3) of the function φ 2 , we have

Hence,

= 0 which has at most two solutions on ]0, S in [ since it is an algebraic equation of degree two in x 1 . Moreover, under (H0)-(H4), the function φ 2 is positive and continuous on [0, S in ] with φ 2 (0) = φ 2 (S in ) = 0. Consequently, φ 2 must have a unique positive extremum on ]0, S in [. Therefore, when (3.6) holds, equation (3.5) may have at most two solutions on [0, S in ] (see Fig. 1(a)).

The following lemma gives the maximal number of extrema of the function F 2 (•) with specific growth functions (2.1).

Lemma Appendix B.2. For specific growth functions (2.1), the function x 1 -→ F 2 (x 1 ) has a unique positive extremum on x 1 1 , x 2 1 .

Proof. From Lemma 3.4, the equation f 2 (S in -x 1 -x 2 , x 1 ) = D defines the function F 2 (•) that can be expressed explicitly by

A straightforward calculation gives the following expression of the derivative of F 2 :