Nico Lingg

Martina Berndtsson

Beate Hintersteiner

Manfred Schuster

Muriel Bardor

Alois Jungbauer
email: alois.jungbauer@boku.ac.at

Highly linear pH gradients for analyzing monoclonal antibody charge heterogeneity in the alkaline range: Validation of the method parameters

Cation exchange chromatography has been routinely used for the quantification of monoclonal antibody (mAb) charge heterogeneity. A previously developed method utilizing pH gradients for the elution instead of salt gradients was validated according to current guidelines proposed by the ICH. The linearity, stability, accuracy, precision and the lower limit of quantification have been determined, using pure charge variant standards. The method is valid for the quantification of mAb samples with a charge heterogeneity between 1% and 50%. Three different approaches to obtaining pure standard material for the validation of bio-analytical methods for the quantification of charge heterogeneity of IgG are presented. These methods are based on salt gradient elution, pH gradient elution and displacement in cation exchange chromatography.

Introduction

Monoclonal antibodies (mAb) are a highly researched class of proteins which represent a very profitable part of the biopharmaceutical industry [START_REF] Elvin | Therapeutic antibodies: market considerations, disease targets and bioprocessing[END_REF]. Indeed, their ability to target a wide range of antigens, e.g. other proteins, carbohydrates, and the strong effector functions they can elicit, e.g. antibody dependent cellular cytotoxicity, complement dependent cytotoxicity, make them an ideal molecule to be used as a therapeutic [START_REF] Schuster | Cancer immunotherapy[END_REF]. Immunoglobulin G (IgG) is the most commonly used isotype. It is a rather large protein with ∼150 kDa and consists of four subunits. In order for the protein to function properly it requires proper folding, correct disulfide bonds and correct N-glycosylation for effector functions [START_REF] Lingg | The sweet tooth of biopharmaceuticals: importance of recombinant protein glycosylation analysis[END_REF][START_REF] Goetze | Assessing monoclonal antibody product quality attribute criticality through clinical studies[END_REF]. Other posttranslational modifications, e.g. methionine oxidation, C-terminal lysine processing, deamidation of asparagine and isomerization of aspartic acid influence the behavior of mAb, as such they can be considered a critical quality attribute (CQA) [START_REF] Liu | Heterogeneity of monoclonal antibodies[END_REF][START_REF] Van Beers | Minimizing immunogenicity of biopharmaceuticals by controlling critical quality attributes of proteins[END_REF]. The emergence of novel purification strategies demands sophisticated analytical methods to ensure comparability of the product [START_REF] Hammerschmidt | Economics of recombinant antibody production processes at various scales: industry-standard compared to continuous precipitation[END_REF][START_REF] Barroso | Functional monolithic platforms: chromatographic tools for antibody purification[END_REF][START_REF] Cramer | Downstream bioprocessing: recent advances and future promise[END_REF][START_REF] Sommer | Combined polyethylene glycol and CaCl2 precipitation for the capture and purification of recombinant antibodies[END_REF][START_REF] Tscheliessnig | Ethanol precipitation for purification of recombinant antibodies[END_REF], and such analytical methods should be validated for their intended purpose [START_REF] Hartmann | Validation of bioanalytical chromatographic methods[END_REF][START_REF] Ich | ICH Harmonised Tripartite Guideline: Validation of Analytical Procedures: Text and Methodology, Q2(R1)[END_REF].

In order to facilitate this goal, we decided to validate our previously presented pH gradient method for the analysis of antibody charge heterogeneity based on cation exchange chromatography (CEX) [START_REF] Lingg | Highly linear pH gradients for analyzing monoclonal antibody charge heterogeneity in the alkaline range[END_REF]. This method is intended to be used for the determination of acidic and basic charge variants, relative to the main charge variant found in a monoclonal antibody batch. In order to validate the method parameters, we purified mAb charge variants, to be used as a standard. The preparative separation of charge variants requires methods with a high resolving power. In this work, we describe different approaches to obtain pure standard material using CEX. The charge variants will be referred to as the main charge variant (MCV) and the acidic charge variants (ACV). The MCV is the main peak eluting in the chromatogram, and is usually not unmodified antibody, it is rather the antibody with the most common modifications [START_REF] Du | Chromatographic analysis of the acidic and basic species of recombinant monoclonal antibodies[END_REF]. The ACV are simply all charge variants eluting before the MCV in CEX. Basic charge variants could not be purified in sufficient amounts to allow usage as standard material, but the method should be valid for all charge variants. Isoelectric focusing (IEF) was used as an orthogonal method to check the purity of the standards.

Our first approach for the preparative scale separation of charge variants was standard CEX elution chromatography, using NaCl as the mobile phase modifier. This method has been used successfully for decades [START_REF] Carta | Protein Chromatography; Process Development and Scale-Up[END_REF] and requires only a simple setup. The next approach was displacement CEX [START_REF] Mcatee | Isolation of monoclonal antibody charge variants by displacement chromatography[END_REF], which has been shown to be a method capable of high resolution separations of mAb charge variants [START_REF] Zhang | Isolation and characterization of therapeutic antibody charge variants using cation exchange displacement chromatography[END_REF]. With the advent of high performance displacers, this approach was very promising. Our last approach focused on the scale up of our own pH gradient elution CEX method into preparative scale, which we have previously shown to have a higher resolution than standard NaCl gradient elution CEX [START_REF] Lingg | Highly linear pH gradients for analyzing monoclonal antibody charge heterogeneity in the alkaline range[END_REF].

In order for this method to be used as a standard for the quantification of mAb charge heterogeneity, we validated the key performance parameters of the method, by following the recommendations for the validation of bioanalytical chromatographic methods presented by Hartmann et al. [START_REF] Hartmann | Validation of bioanalytical chromatographic methods[END_REF], while using the terminology from Peters et al. in this work [START_REF] Peters | Validation of new methods[END_REF]. In detail the following parameters have been determined: selectivity, calibration model (linearity), stability, accuracy (bias), precision (repeatability, intermediate precision) and the lower limit of quantification (LLOQ). This was done in a similar fashion to the work published by Tscheliessnig and Jungbauer [START_REF] Tscheliessnig | High-performance monolith affinity chromatography for fast quantitation of immunoglobulin G[END_REF] and Rozhkova [START_REF] Rozhkova | Quantitative analysis of monoclonal antibodies by cationexchange chromatofocusing[END_REF]. The upper limit of quantification (ULOQ) was not determined, instead a charge heterogeneity of 50% was taken as the upper limit.

Our final goal was to provide evidence that our previously presented method is indeed valid for the quantification of mAb charge variants.

Other methods, such as imaging capillary isoelectric focusing (iCIEF), can also be used for quantification of mAb charge variant content and are well established in the biopharmaceutical industry [START_REF] Anderson | Applications of imaged capillary isoelectric focussing technique in development of biopharmaceutical glycoprotein-based products[END_REF][START_REF] Sosic | Application of imaging capillary IEF for characterization and quantitative analysis of recombinant protein charge heterogeneity[END_REF]. Such methods have shown great inter-laboratory robustness and are reliable methods [START_REF] Salas-Solano | Robustness of iCIEF methodology for the analysis of monoclonal antibodies: an interlaboratory study[END_REF].

Material and methods

Materials and chemicals

All chemicals were of analytical grade, unless stated otherwise. 3-Morpholino-2-hydroxypropanesulfonic acid (MOPSO), 4-(2hydroxyethyl)piperazine-1-ethanesulfonic acid (HEPES), N,Nbis(2-hydroxyethyl)glycine (bicine), 3-(cyclohexylamino)-2-hydroxy-1-propanesulfonic acid (CAPSO) and 3-(cyclohexylamino)-1-propanesulfonic acid (CAPS), 2-(N-Morpholino)ethanesulfonic acid (MES), NaH 2 PO 4 , Tergitol (70% in water) and isopropanol were purchased from Sigma-Aldrich (St. Louis, MO, USA). NaCl,

Preparative chromatography

The mAb used for the validation was a chimeric IgG 1, provided by Apeiron Biologics (Vienna, Austria). An ÄKTA Explorer 100 (GE Healthcare) was used for all preparative chromatography steps. The outlet was monitored at 280 nm.

Salt gradient elution

A Dionex ProPac WCX-10 column, 9 × 250 mm (Thermo Fisher, Sunnyvale, CA, USA) was used. The flow rate was 2.5 ml/min (236 cm/h). The buffers used for the salt gradient experiments are described in Table 1. The method used consisted of a 1 CV equilibration step at 0% B, 10 ml injections via sample loop, a 1 CV wash step at 0% B, a 1 CV gradient to the initial elution conductivity and the elution gradient with varying gradient length, found in Table 1.

Displacement chromatography

For displacement chromatography, a Mono S column, 4 × 200 mm (GE Healthcare, Uppsala, Sweden) was used. The flow rate was 0.25 ml/min (76 cm/h). Before each experiment, the dynamic binding capacity (DBC) at 10% breakthrough was determined for each mobile phase. The buffers and conditions used for displacement chromatography can be found in Table 2. The displacement experiments were performed according to Zhang et al. [START_REF] Zhang | Isolation and characterization of therapeutic antibody charge variants using cation exchange displacement chromatography[END_REF] and McAtee and Hornbuckle [START_REF] Mcatee | Isolation of monoclonal antibody charge variants by displacement chromatography[END_REF].

pH gradient elution

For salt gradient and pH gradient preparative experiments, a Dionex ProPac WCX-10 column, 9 × 250 mm (Thermo Fisher, Sunnyvale, CA, USA) was used. The flow rate for the ProPac column was 2.5 ml/min (236 cm/h). The buffers used for pH gradient experiments are found in Table 3. The various gradients used are detailed in Table 4. The collected acidic and main fractions from each run were concentrated and buffer exchanged with Amicon Ultra-15 Centrifugal Filter Units, 50 kDa (Merck Millipore, Darmstadt, Germany). The diafiltration buffer was 20 mM histidine, 0.1% Tween 20, pH 6.0. The gradient for the second purification round of the ACV standard consisted of two segments, a steeper initial gradient and a shallower gradient to separate the ACV from the main charge variant. The concentrated and buffer exchanged mAb was adjusted to a concentration of 2 mg/ml and used as the MCV and ACV standards.

Isoelectric focusing

Immobilized pH gradient (IPG)-polyacrylamide-gels (size: 125 × 260 × 1 mm) were cast on a GELbond-PAG-film backing (GE Healthcare), polymerized at 50 • C for 2 h, washed and dried following a procedure previously described by Westermaier [START_REF] Westermeier | Electrophoresis in Practice: A Guide to Methods and Applications of DNA and Protein Separations[END_REF]. The desired pH gradient from 7.0 to 11.0 was obtained by graphic interpolation from a recipe previously published by Gorg et al. [START_REF] Gorg | Very alkaline immobilized pH gradients for two-dimensional electrophoresis of ribosomal and nuclear proteins[END_REF].

Directly before use, the gel was cut in 3 mm strips, which were rehydrated for 2 h in a solution containing 6 M urea, 2% tergitol, 2% pharmalyte 3-10, 10% glycerol and 16% isopropanol.

The rehydrated strips were put onto the IPGphor (GE Healthcare) using the manifold. The electrodes were positioned on the acidic and basic ends of the strips with wetted paper wicks between gel and electrode to ensure good contact. The samples, 15 g of mAb, were applied by cup loading under a covering layer of paraffin and electrophoresis was performed overnight (at 150 V for 1 h, followed by 300 V for 3 h, 3500 V for 18 h). After electrophoresis, the strips were stained with Coomassie blue G250.

Analytical chromatography

An Agilent 1220 Infinity LC System (Santa Clara, CA, USA) equipped with a SIM sample-cooler (Scientific Instruments Manufacturer GmbH, Oberhausen, Germany) was used for all validation experiments. The column used was a Dionex ProPac WCX-10, 4 × 250 (Thermo Fisher). The column oven and the auto-sampler temperatures were set to 30 • C and 10 • C, respectively. These temperatures were chosen due to instrument limitations.

The samples for the calibration model, the accuracy, precision and the stability experiments were prepared by mixing appropriate volumes of ACV and MCV standard. The outlet was monitored at 280 nm. The flow rate was 1 ml/min (477.5 cm/h). The buffers were the same as those used for the preparative experiments in Table 3. A gradient from 35% to 55% buffer B in 8 CV was used. The chromatograms were then integrated with the Chemstation software suite (Agilent).

Results and discussion

Creation of charge variant standard

In order to validate the method parameters, pure IgG charge variant standards were necessary. For such preparations, different CEX approaches were conceived and tested. The first approach was cation exchange chromatography combined with salt gradient elution.

Salt gradient elution

Both the Mono S and the ProPac WCX-10 columns were tested, but the latter column was used for further optimization due to the higher resolution, owing to the faster mass transfer.

The method was optimized by testing elution behavior at different pH values and adjusting the salt gradient, the conditions tested for the optimization can be found in Table 1. The elution was optimized by running increasingly narrow gradients around the elution peaks of interest, thereby increasing the resolution. Optimization experiments were performed with 190 g/ml injections, with a total volume of 10 ml, but subsequent scale up to 1.9 mg/ml was unsuccessful. Fig. 1 shows the optimized salt gradient with 190 g/ml (Fig. 1A) and 950 g/ml (Fig. 1B) injections. A drastic decrease in resolution was observed even under relatively low loading amounts (<1 mg/ml stationary phase). Due to the low binding capacity observed at pH 7, we did not attempt separation at a higher pH. Since scale up of salt gradient elution failed, we pursued displacement chromatography as a method for creating the standard.

Displacement chromatography

Displacement chromatography was tested at three different conditions, as described in Table 2. Expell SP1 was chosen as the displacer, as it has proven to be useful in IgG separations [START_REF] Mcatee | Isolation of monoclonal antibody charge variants by displacement chromatography[END_REF][START_REF] Zhang | Isolation and characterization of therapeutic antibody charge variants using cation exchange displacement chromatography[END_REF][START_REF] Khawli | Charge variants in IgG1: isolation, characterization, in vitro binding properties and pharmacokinetics in rats[END_REF]. No separation of charge variants was detectable after analysis with IEF. A representative result from the pH 7.0 displacement run is shown in Fig. 2. Due to the disappointing results in all three displacement experiments, we decided to scale up our pH gradient elution method to a semi-preparative scale.

pH gradient elution

Finally preparative separation of the charge variants using pH gradient elution was used, lacking an orthogonal method for creating the standard. The initial mAb sample was purified as shown in Table 4 and acidic and main peak fractions were collected, as shown in Fig. 3(A). After concentration and buffer exchange of the eluate with ultra-/diafiltration, the standards were used in a second round of purification using a shallower gradient to remove the remaining unwanted charge variants. All three chromatograms can be seen in Fig. 3, including the fractionation that was used to obtain the pure standard material. It should be noted that the linearity of the pH gradients stays almost unaffected, even with higher loads of antibody. The gradients used can be found in Table 4.

The standards obtained this way will be referred to as main charge variant (MCV) and acidic charge variants (ACV). Because of non-orthogonality of the method used for creating the standard, IEF was used to confirm that the standard material indeed consisted of purified and separated charge variants. Fig. 4 shows the original mAb sample and both of the standards used in this work. No main band can be detected in the ACV sample, and only a faint band with a lower isoelectric point (pI) is visible in the MCV sample. It was concluded that these two samples will be of sufficient purity to validate the analytical method.

Testing of the standard material

The purified MCV and ACV standards were analyzed using the analytical pH gradient elution method to check if the standards are of sufficient purity to be used for the validation and compared with the original mAb sample. A gradient from 35% to 55% buffer B (pH ∼8.9-9.4, see Table 3 for composition) in 8 CV was used. This specific gradient was able to separate the charge variants, as shown in Fig. 5, while allowing to reduce the run time. Therefore, it was used for the whole validation. Both standards are of over 99% purity and elute at different times, as shown in Fig. 5C. It was concluded from these experiments that MCV and ACV will be suitable standard material for the subsequent validation experiments.

Calibration model

In order to confirm the linearity of the calibration model and to find the linear range of the method, spiking experiments were performed. MCV was spiked with ACV, to obtain relative abundances from 1% to 50% acidic species. The calibration model was tested by fitting a linear model to the data obtained with eight concentration levels evenly spaced across the calibration range with three replicates [START_REF] Peters | Validation of new methods[END_REF]. The results can be found in Table 5 and Fig. 6. The results of the linear regression gave a slope of 157.9 and an intercept of 143.0 with a correlation of regression of 0.9996. The residuals were equally (Fig. 7) distributed and the RSD was lower than 10% over the whole range. One value with an unusually large residual was observed, but Grubb's outlier test (˛= 5%) determined that this was not an outlier [START_REF] Hartmann | Validation of bioanalytical chromatographic methods[END_REF]. ANOVA was used to confirm the validity of a linear calibration model, with the following results: linearity F test (˛= 1%), F* = 16640.41 > F = 7.95, which leads to the acceptance of linear model; lack of fit F test (˛= 1%), F* = 1.23 < F = 4.20, which leads to the acceptance of the null hypothesis of good model fit.

Accuracy and precision

The accuracy and bias were determined by the statistical evaluation of duplicate experiments performed on eight days at three different concentrations, as recommended by Peters et al. and Hartmann et al. [12,[START_REF] Peters | Validation of new methods[END_REF]. The concentration levels selected were the LLOQ, the ULOQ and the median, i.e. 1%, 50% and 25% ACV, respectively. The precision was determined as the overall RSD, with acceptance criteria being 20% RSD for the lower limit of quantification and 15% for the other two spiking levels. The accuracy was determined as the deviation of the grand means from the respective reference values, with acceptance criteria being ±20% bias for the lower limit of quantification and ±15% for the other two spiking levels [START_REF] Peters | Validation of new methods[END_REF][START_REF] Shah | Bioanalytical method validation-a revisit with a decade of progress[END_REF][START_REF] Shah | Analytical methods validation: bioavailability, bioequivalence and pharmacokinetic studies[END_REF]. The results from measuring at the LLOQ, the median and the ULOQ can be seen in Table 6. Both the bias and the RSD are within the acceptance criteria. The bias at the LLOQ is rather high, with -17.6%, suggesting that quantification of charge variants below 1% in relative abundance will not be practical. This also served as confirmation of our estimation of the LLOQ for this method.

Selectivity

We assume that this method will be used for the determination of charge heterogeneity of purified monoclonal antibodies only, and as such the components that might be expected to be present in the samples analyzed are small molecules, such as buffering compounds, amino acids, sugars and detergents. As such, selectivity, "the ability to measure the analyte in the presence of components which might be expected" [START_REF] Peters | Validation of new methods[END_REF], has been shown previously [START_REF] Lingg | Highly linear pH gradients for analyzing monoclonal antibody charge heterogeneity in the alkaline range[END_REF].

Stability

The stability of mAb samples inside the autosampler, at 10 • C, and the stability of the pH gradient buffers was tested over a time frame of 48 h. Samples were measured in triplicates at 0, 12, 24 and 48 h. The distribution of the measured values over time can be seen in Fig. 8. Those results demonstrated that the samples and the method are stable for 48 h and the method can therefore be used over night or over the weekend.

Conclusions

We presented the development of a highly linear pH gradient elution method in our previous work [START_REF] Lingg | Highly linear pH gradients for analyzing monoclonal antibody charge heterogeneity in the alkaline range[END_REF]. In the present complementary work, we have validated key method parameters, such as the accuracy and precision. The method uses commercially available stationary and mobile phase components, which allows others to implement this method easily. Other methods such as IEF provide a higher resolution than our method, but can only be considered as semi-quantitative method, cannot be automated as easily and cannot be scaled up to semi-preparative scale. The method for the quantification of mAb charge variants is valid for the quantification of purified mAb samples with a charge variant content as low as 1%. The linearity of the calibration range was confirmed through statistical testing and also showed an excellent coefficient of correlation of 0.9996. 1% charge variants was confirmed as the LLOQ, with an accuracy (bias) of -17.6% and a precision of 8.9%. For samples measured at the median or the ULOQ, the method had excellent accuracy and precision of below 2%. We have previously presented methods for the creation of gradients from pH 7.0 to pH 10.5. The big advantage of our method is that the linearity of the gradient is still high even when a narrower gradient is being used to increase the resolution. This means that our method can be easily customized for the analysis of different mAb, just by varying the initial and the final buffer concentration. In this study, we have used an antibody with a chromatographic pI of 9.2 and have adjusted our method gradient from 35% to 55% buffer B accordingly, which corresponded to a pH gradient from 8.9 to 9.4.

Compared to iCIEF, our method shows comparable precision [START_REF] Sosic | Application of imaging capillary IEF for characterization and quantitative analysis of recombinant protein charge heterogeneity[END_REF], but no thorough validation could be found in which accuracy of iCIEF was determined. We conclude that our method is comparable to other standard methods. The biggest difficulty in performing this study was the production of the standard material. Our attempts at using standard salt gradient elution cation exchange chromatography proved unfruitful. We were unable to find conditions under which we had a sufficient resolution with an adequate protein load. The approach presented by Rozhkova [START_REF] Rozhkova | Quantitative analysis of monoclonal antibodies by cationexchange chromatofocusing[END_REF], i.e. cation exchange chromatography combined with salt gradient elution, only worked for very small loading amounts for our mAb and additionally did not achieve the resolution we observed in pH gradient elution experiments. Unfortunately Rozhkova does not present analytical data of her standard material as well as more detailed information about loading amounts. Accordingly it is difficult to compare the two methods.

Our next approach, displacement chromatography [START_REF] Mcatee | Isolation of monoclonal antibody charge variants by displacement chromatography[END_REF], did not yield a separation of charge variants, under any conditions. This is in stark contrast to the results that Zhang et al. [START_REF] Zhang | Isolation and characterization of therapeutic antibody charge variants using cation exchange displacement chromatography[END_REF] have presented previously. Eventually, we used a preparative pH gradient elution method to prepare our standards, and confirmed their identity as charge variants by IEF. This approach proved to be highly effective in separating large amounts. It has previously been shown that induced pH gradients can be used for the preparative separation of mAb charge variants [START_REF] Vetter | Mixed-beds of strong and weak anion exchange resins for protein separations with step-induced pH gradients[END_REF][START_REF] Vetter | Predicting retention and resolution of protein charge variants in mixed-beds of strong and weak anion exchange resins with step-induced pH gradients[END_REF][START_REF] Pabst | Separation of protein charge variants with induced pH gradients using anion exchange chromatographic columns[END_REF], and we are happy to report that our external pH gradients are able to achieve the same feat. Solubility problems due to elution close to the pI were not observed, even at the preparative scale used in this work, making this a high resolution method for the preparative separation of protein charge variants.

Fig. 1 .

 1 Fig.1. Chromatogram of the semi-preparative separation of charge variants using salt gradient elution. Panel (A) shows a load of 3 mg (∼190 g/ml column) and panel (B) shows a load of 15 mg (∼950 g/ml column). A drastic decrease in resolution is apparent at higher loading. Stationary phase: Dionex ProPac WCX-10 column, 9 × 250 mm. Flow rate: 2.5 ml/min (236 cm/h). Mobile phase: 10 mM phosphate, pH 7.0, buffer A with 50 mM NaCl, buffer B with 100 mM NaCl. Gradient from 60 to 65.5 mM NaCl in 5 CV.

Fig. 2 .

 2 Fig. 2. IPG IEF of the displacement fractions from a displacement experiment on a Mono S column, 5 × 200 mm. Flow rate: 0.25 ml/min (76 cm/h). Mobile phase: 30 mM MOPSO, 15 mM NaCl, pH 7.0. Displacer: 5 mM Expell SP1. Lanes 1-10 contain displacement fractions with higher numbers denoting later eluting fractions. O is the initial sample. No difference in any fraction could be observed. No marker was used, since only the isoelectric point (pI) relative to the initial sample was relevant.

Fig. 3 .Fig. 4 .

 34 Fig.3. Chromatograms of the preparative pH gradient elution separation of charge variants on a ProPac WCX-10, 9 × 250 mm column. Flow rate: 2.5 ml/min (236 cm/h). Mobile phase A: 5.5 mM HEPES, 4.2 mM Bicine, 9.5 mM CAPSO, 0.8 mM CAPS, 6.3 mM NaCl, pH 8.0, mobile phase B: 10.5 mM Bicine, 2.5 mM CAPSO, 7.0 mM CAPS, pH 10.5. Panel (A) shows the initial purification and fractionation for the acidic (blue) and the main (orange) standard. Graqdient from pH 8.0 to 9.8 in 14 CV. After concentration and buffer exchange those fractions were used in a second round of purification to further increase the purity. Panels (B) and (C) show the purification and fractionation of the second round for acidic and main standard, respectively. For the acidic standard a 2 CV gradient from pH 8.0 to 8.9 was followed by a 10 CV gradient to pH 9.4. For the main standard a 10 CV gradient from pH 8.9 to 9.4 was used. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 5 .Fig. 6 .

 56 Fig. 5. Evaluation experiments for the suitability of the standard. The MCV (orange) and ACV (blue) standards were tested using the pH gradient elution method. Stationary phase: Dionex ProPac WCX-10, 4 × 250 mm. Flow rate: 1 ml/min (477.5 cm/h). Mobile phase A: 5.5 mM HEPES, 4.2 mM Bicine, 9.5 mM CAPSO, 0.8 mM CAPS, 6.3 mM NaCl, pH 8.0, mobile phase B: 10.5 mM Bicine, 2.5 mM CAPSO, 7.0 mM CAPS, pH 10.5. A gradient from 35% to 55% B in 8 CV was used. Panel (A) and (B) show the chromatograms of the MCV and ACV standard, respectively. The gray dashed line represents the pH at the column inlet, not a measured signal. Panel (C) shows the overlay of both chromatograms from 5 to 25 min. The samples are of high purity and are suited for the use as a charge variant standard. The small initial peak is a solvent peak, consisting of the histidine and Tween 20 of the mAb buffer. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 8 .

 8 Fig. 8. Results from the stability experiments over 48 h. Samples with 1%, 25% and 50% ACV were measured in triplicates at 0, 12, 24 and 48 h.

Table 1

 1 Summary of the salt gradient elution experiments performed on a semi-preparative ProPac WCX-10 column. [M] is the modifier (NaCl) concentration, with the indices i and f meaning initial and final respectively.

	Buffer	[NaCl] buffer A (mM)	[NaCl] buffer B (mM)	pH	[M]i (mM)	[M]f (mM)	Gradient slope (CV)
	30 mM MES	0	1000	6.0	50	350	20
	10 mM phosphate	50	100	7.0	50	100	5
	10 mM phosphate	50	100	7.0	70	77	5
	10 mM phosphate	50	100	7.0	63	70	5
	10 mM phosphate	50	100	7.0	63	65.5	5
	10 mM phosphate	50	100	7.0	60	65.5	5

Table 2

 2 Summary of the displacement chromatography experiments performed on a Mono S column.

	Buffer	pH	Displacer (mM)	Load (mg/ml)	Load (% DBC10)
	40 mM MES, 15 mM NaCl	5.5	5	35	50
	30 mM MES, 100 mM NaCl	6.1	5	14.9	60
	30 mM MOPSO, 15 mM NaCl	7.0	5	33.6	60

Table 3

 3 The buffer system used throughout this work for the creation of linear pH gradients.

		HEPES	Bicine	CAPSO	CAPS	NaCl	pH
	Buffer A [mM]	5.5	4.2	9.5	0.8	6.3	8.0
	Buffer B [mM]	0.0	10.5	2.5	7.0	0.0	10.5
	Na 2 HPO 4 , NaOH, Urea and glycerol was purchased from Merck
	(Darmstadt, Germany). Pharmalyte 3-10 was purchased from (GE
	Healthcare, Uppsala, Sweden). Histidine, research grade (>98.5%
	purity) was from Serva (Heidelberg, Germany). Tween 20 was
	obtained from Bio-Rad (Hercules, CA, USA). Expell TM SP1 (Sachem
	Inc., Austin, TX, USA) was used as the displacer in displacement
	chromatography experiments.

Table 4

 4 Outline of the pH gradient elution method used for the purification of the standard material.

	Experiment	Wash (CV)	Initial pH	Gradient (CV)	Final pH	Gradient (CV)	Final pH
	1st round	1	8.0	14	9.8
	2nd round, acidic	1	8.0	2	8.9	10	9.4
	2nd round, main	1	8.9	10	9.4

Table 5

 5 Calibration results.

	ACV content (%)	Response (AU*s)	RSD (%)
	0	97.2	7.7
	1	280.2	6.7
	3	632.2	1.4
	5	945.3	0.3
	10	1719.7	0.6
	20	3406.9	0.6
	40	6355.5	4.4
	50	8079.8	1.0

Table 6

 6 Results of the accuracy and precision measurement. Fig. 7. Residuals of the calibration curve.

	Relative amount of	Average relative charge variant	Bias (%)	Standard deviation	Relative standard deviation (%)
	spiked ACV (%)	content measured (%)
	1	0.82	-17.6	0.07	8.9
	25	25.37	1.5	0.38	1.5
	50	49.96	-0.1	0.32	0.6

Acknowledgments

We would like to thank Dr. Anne Tscheliessnig for her help with the statistical analysis and Nikolaus Hammerschmidt and Peter Satzer for their critical review of the manuscript. This work was supported by the Austrian Research Promotion Agency (FFG). Nico Lingg was a fellow from the A*STAR Research Attachment Programme (ARAP). The collaboration between the Glyco-MEV, University of Rouen, France and BTI, A*STAR, Singapore was supported by the Merlion 2011 initiative program called Glyco-TOOLS. We would like to thank GE Healthcare for providing the Mono S columns used in this study.


T Guérin

M Dolgushev

O Bénichou

R Voituriez

Imperfect Narrow Escape problem

We consider the kinetics of the imperfect narrow escape problem, i.e. the time it takes for a particle diffusing in a confined medium of generic shape to reach and to be adsorbed by a small, imperfectly reactive patch embedded in the boundary of the domain, in two or three dimensions. Imperfect reactivity is modeled by an intrinsic surface reactivity κ of the patch, giving rise to Robin boundary conditions. We present a formalism to calculate the exact asymptotics of the mean reaction time in the limit of large volume of the confining domain. We obtain exact explicit results in the two limits of large and small reactivities of the reactive patch, and a semi-analytical expression in the general case. Our approach reveals an anomalous scaling of the mean reaction time as the inverse square root of the reactivity in the large reactivity limit, valid for an initial position near the extremity of the reactive patch. We compare our exact results with those obtained within the "constant flux approximation"; we show that this approximation turns out to give exactly the next-to-leading order term of the small reactivity limit, and provides a good approximation of the reaction time far from the reactive patch for all reactivities, but not in the vicinity of the boundary of the reactive patch due to the above mentioned anomalous scaling. These results thus provide a general framework to quantify the mean reaction times for the imperfect narrow escape problem.

I. INTRODUCTION

How much time does it take for a random walker to reach a target point? The answer to this question has received a lot of attention in the last decade in the physics literature [START_REF] Redner | A guide to First-Passage Processes[END_REF][START_REF] Condamin | [END_REF][3][4][5][6][START_REF] Metzler | First-passage phenomena and their applications[END_REF][START_REF] Schuss | [END_REF][9]. First passage problems appear in various areas of biological and soft matter physics and are in particular relevant to the problem of reaction kinetics, since two reactants have to meet before being able to react [START_REF] Rice | Diffusion-Limited Reactions[END_REF][START_REF] Berg | [END_REF]. When the reaction is "perfect", i.e. when it occurs instantaneously upon each encounter, its kinetics is controlled by the first passage statistics of one reactant molecule, seen as a random walker, to the second reactant, seen as a "target". However, many reactions do not occur at first contact between the random walker and the targets, leading to imperfect reactivity. Imperfect reactivity can have diverse origins at the microscopic scale, such as orientational constraints on the reactive particles [START_REF] Berg | [END_REF], the fact that the surface of the reactive particles is not entirely covered by reactive patches (such as in the chemoreception problem [12]), the need to overcome an energetic [13] (or entropic [14]) activation barrier before reaction, the presence of a gate that can be randomly closed or opened when the reactant meets the target [15,16], etc (see Ref. [START_REF] Grebenkov | Chemical Kinetics: Beyond the Text[END_REF] and references therein for a recent review on imperfect reactivity).

Imperfect reactivity was early investigated for molecules diffusing in infinite space [START_REF] Berg | [END_REF][START_REF] Collins | [END_REF][19][20] (with an imposed concentration at infinity). The search problem for a single random walker moving in a confining volume for an imperfect target, initially considered in Ref. [21] for centered spheres, has also attracted recent attention and several asymptotic results for imperfect search kinetics have been derived [22][23][24][25][26][27][28][29][30]. Recently, explicit asymp-totics of the reaction time statistics have been obtained for general Markovian random walks [31]. Besides the case of reactive targets located in the bulk of a confining domain, the narrow escape problem (NEP) consists in calculating the escape time of a random walker out of a confining domain, through a small window at the boundary of the domain (see Figure 1(a)). While the NEP is now well characterized for perfect reactions, for spherical domains [START_REF] Schuss | [END_REF][32][33][34][35] and large domains of arbitrary shapes [36], fewer results are available for imperfect reactions (i.e. for a partially adsorbing patch). The imperfect narrow escape problem has been investigated for particular geometries in cylindrical [37,38] or spherical domains [39] [40] in which case the analysis depends on the eigenfunctions of the particular confining volume that is considered and relies on the so-called uniform flux approximation introduced in Ref. [20].

The aim of the present paper is to apply the formalism introduced in Ref. [31] to cover the case of the imperfect narrow escape problem in a domain of generic shape. Our formalism is asymptotically exact in the limit of large confining volume -it does not involve the constant flux approximation -and provides explicit results in both regimes of small and large reactivity. Of particular interest for imperfect reaction problems is the mean reaction time when the initial position is located on the reactive patch; this time is exactly zero for perfect reactions and scales as 1/κ for targets in the interior of the volume. We identify a region for which the reaction time behaves anomalously with the reactivity κ. This region, which does not exist for targets in the bulk of the confining domain, is located at the boundary of the imperfectly reactive patch. While one would naively expect this time to be inversely proportional to κ, we find in- stead that when the initial position is at the boundary of the reactive domain, the mean reaction time T e is actually anomalously high (∝ κ -1/2 ) and follows the exact asymptotics

T e ∼ κ→∞ V (d -1)(2π κDa) 1/2 a d-2 .
(

Here, d = 2 or d = 3 is the spatial dimension, D is the diffusion coefficient, a is the radius of the reactive patch, and V the volume of the confining domain. Here, we assume that the confining volume is taken large enough, and the patch small enough, so that the confining boundaries at the vicinity of the target can be considered as a flat wall in which the reactive patch is embedded, the latter being considered as a line segment of length 2a in d = 2 or as a flat disk of radius a in d = 3. We show below how to obtain this anomalous scaling relation by solving a Wiener-Hopf integral equation. We will also show how this "anomalous" behavior (1) of T with κ can be related to the divergence of fields in Laplacian problems near surfaces presenting asperities, as occurs in electrostatics near conducting edges [START_REF] Jackson | Classical electrodynamics[END_REF] or in the coffee ring effect [START_REF] Deegan | [END_REF].

More generally, we show that the mean escape time for an arbitrary initial position far from the target, and for any finite reactivity κ satisfies the following exact asymptotics:

T /V ∼ r a 1 πD ln(r/a) + C ∞ (d = 2) -1 2πDr + C ∞ (d = 3) , (2)
where C ∞ is independent of the initial distance r from the target. For finite values of the reactivity κ, we show that C ∞ can be obtained through a semi analytical procedure.

In the limit of large reactivity, we show that C ∞ can be determined explicitly and is given by :

C ∞ ∼ κa D ln 2 πD + 1 π 2 κa ln 8κa D + γ e + 1 (d = 2) 1 4Da + 1 4πκa 2 ln 2κa D + γ e + 1 (d = 3) (3)
where γ e is Euler's constant. This expression is understood as the first two terms in the expansion of C ∞ in powers of 1/κ. Interestingly, this result shows that the term C ∞ is not analytic in powers of κ, which originates from the anomalous scaling [START_REF] Redner | A guide to First-Passage Processes[END_REF]. Finally, we also give exact results in the small reactivity limit, which will be found to be exactly the same (at first order) as the results obtained within the self-consistent, "constant flux approximation" that has been invoked to study the imperfect narrow escape in the literature [20,37,39,40]. It is found that, far from the reactive patch this approximation is very accurate (for any reactivity), while it fails for initial positions close to the reactive patch, and in particular does not predict the behavior (1) near the boundary of the patch. The outline of the paper is as follows. First, we recall the formalism of Ref. [31] in the particular case of the imperfect narrow escape problem for diffusing particles to obtain equations for the mean escape time in the large volume limit (Section II A). We show how the formalism can be presented under the form of an integral equation that is suitable for studying the large and small reactivity limits in Section II B. The large reactivity limit is investigated in Section III where Eqs. ( 1) and (3) are derived. In this Section, we also give a simple scaling argument that relates the anomalous behavior of T near the extremity of the patch to the divergence of electric fields near the edges of conducting objects. The small reactivity limit is examined in Section IV. Last, we study briefly how the constant flux approximation can be implemented within our formalism in Section V. An exact, but formal solution for any reactivity parameter (that requires numerical tools, however) is presented in Appendix A.

II. FORMALISM FOR THE IMPERFECT NARROW ESCAPE PROBLEM IN THE LARGE VOLUME LIMIT

A. General formalism

We consider the stochastic motion of an overdamped particle moving with diffusion coefficient D in a confining volume Ω. The boundary of the volume is ∂Ω and contains a small window S r which is partially reactive, the rest of the confining boundary is assumed to be smooth and perfectly reflecting, see Fig. 1(a). We assume that S r is formed by the region of the surface at geodesic distance less than a from the center, and a is called the radius of the patch S r . The Fokker-Planck equation for the probability density p(r, t) to observe the particle at position r and time t is

∂ t p = D ∇ 2 p (r ∈ Ω), (4) n • ∇p = 0 (r ∈ ∂Ω\S r ), (5) D n • ∇p + κ p = 0 (r ∈ S r ). ( 6
)
where n is the unit vector normal to the surface, pointing to the exterior of the volume. For a partially reactive surface, the reactivity parameter κ is defined in such a way that the probability that the particle is absorbed by an infinitesimal surface element dS located around r s during dt is κ p(r s , t)dSdt. We assume that the space dimension is d = 2 (2D) or d = 3 (3D). It is very well known that an equation for the mean first passage time can be obtained by identifying the adjoint transport operator [START_REF] Redner | A guide to First-Passage Processes[END_REF], which in our case leads to the following equation for the mean reaction time T (r) to the target, where r now represents the initial position of the particle:

D∇ 2 T = -1 (r ∈ Ω), ( 7
) n • ∇ T = 0 (r ∈ ∂Ω\S r ), ( 8
) D n • ∇ T + κ T = 0 (r ∈ S r ). (9)
Integrating Eq. ( 7) over the whole volume, and using the divergence formula and the boundary conditions leads to the exact integral relation:

κ Sr dS(r) T = V, (10)
where V = |Ω| is the volume of the domain. Now, we consider the large volume limit, which is obtained when the confining volume extends without changing its shape, keeping constant the size of the target and the initial distance to the target. We define the rescaled mean escape time Φ by

Φ(r) = lim V →∞ T (r) /V. (11)
In the large volume limit, the boundary at the vicinity of the reactive target becomes increasingly similar to a flat surface in which the reactive patch is a flat disk of radius a in 3D (or a flat segment of length 2a in 2D). Here, we denote the distance to the reflecting surface containing reactive patch as z, see Fig. 1(b). With this in mind, inserting the ansatz (11) into the above equations yields a closed system of equations in the large volume limit:

∇ 2 Φ = 0 (if z > 0), (12)
κ Sr dS(r)Φ = 1, ( 13
)
D∂ z Φ = 0 (if z = 0, |r| > a), κ Φ (if z = 0, |r| < a). (14)
Importantly, we see that in the large volume limit, Eq. ( 11), the obtained equations are independent of the shape of the confining volume, which is present only though the scale factor V in the definition of Φ. We have directly controlled this aspect by performing numerical stochastic simulations of trajectories in the confined domain. The results of such simulations are shown on Fig. 2 and confirm that our formalism correctly predicts the mean reaction time in the large volume limit, independently on the shape of the confining domain. Equations ( 12), ( 13), ( 14) generalize the formalism of Ref. [36] to the case of imperfect reactions. The fact that the mean first reaction time scales with the volume is actually more general than the specific diffusive walk that we have considered here [31]. To solve the above equations, we may be tempted to use spheroidal coordinates, which can be used to solve the problem for either infinite or vanishing reactivity. For finite reactivity however, the resulting equations in such coordinates involve Robin conditions with non-uniform coefficients, so that the mean reaction time can be obtained only in terms of the solution of an infinite linear system. This procedure is described in Appendix A, and it indeed leads to a generic numerical solution that will be useful to test our analytical insights in all the paper. However, it is not suitable for analytical calculations. For this reason, we adopt a different approach, consisting in deriving an integral equation satisfied by Φ on the reactive patch.

B. Obtaining an integral equation for the mean reaction time

Let us first characterize the large distance behavior of the rescaled mean first reaction time Φ. The condition (13), combined with the boundary condition at z = 0, implies that

S0 dS n • ∇Φ = 1/D, (15)
for any surface S 0 whose intercept with the plane z = 0 encloses the reactive patch. Taking such surface S 0 to be a half-disk of radius R (in 2D) or a half-sphere (in 3D), we see that

∂ r Φ ∼ r→∞ 1/[πDr] (d = 2), 1/[2πDr 2 ] (d = 3), ( 16
)
where r is the distance to the center of the reactive patch. Hence, the behavior of Φ for large r takes the form

Φ(r) ∼ |r|→∞ 1 πD ln |r| + C ∞ + o(1) (d = 2), -1 2πDr + C ∞ + o(1) (d = 3), ( 17
)
where C ∞ does not depend on r. The quantity C ∞ thus characterizes the behavior of the mean reaction time far from the target, it could be used in matched asymptotics expansions if one aims to identify the first passage times distributions, as in ref. [22] (which deals with interior targets). (parameters are indicated in the legend) compared to general theoretical expressions as obtained in Appendix A. In all simulations, we used a time step ∆t = 10 -4 a 2 /D. Boundary conditions are implemented as follows: if, at the end of a time step, the random walker falls outside the domain, then if it is "behind" a reflecting wall it is reflected with respect to this wall, and if it falls "behind" the reactive patch, it is absorbed with probability Pa = κ πdt/D (in which case the trajectory ends) and it is reflected with probability 1 -Pa, see Ref. [43].

Let us now introduce the Green's function G N for the Laplace problem with Neumann boundary conditions at z = 0 (including on the reactive region). Such Green's function satisfy

∇ 2 G N (r|r 0 ) = -δ(r -r 0 ), (18) n • ∇G N (r|r 0 ) = 0 (z = 0). ( 19
)
The expression of G N is easily found by using the image method [START_REF] Barton | Elements of Green's functions and propagation[END_REF]:

G N (r|r 0 ) = -1 2π [ln |r -r 0 | + ln |r -r * 0 |] (d = 2), 1 4π 1 |r-r0| + 1 |r-r * 0 | (d = 3), (20)
where r * 0 represents the symmetric image of r 0 with respect to the plane z = 0. We note that the large r be-

havior of G N is G N (r|r 0 ) = |r|→∞ -1 π ln |r| + o(1/r) (d = 2), 1 2πr + o(1/r 2 ) (d = 3), (21)
We now use manipulations that are standard in Green's function problems [START_REF] Barton | Elements of Green's functions and propagation[END_REF] to put the problem for Φ on the form of an integral equation. Using the Eq. ( 18) and ∇ 2 Φ = 0 we see that the following equality holds:

Φ(r 0 ) = z≥0 dr[G N (r|r 0 )∇ 2 r Φ -Φ(r)∇ 2 r G N (r|r 0 )]. (22)
Using the divergence formula, we obtain

Φ(r 0 ) = S dS(r)n • [G N (r|r 0 )∇ r Φ(r) -Φ(r)∇ r G N (r|r 0 )], ( 23
)
where S is any closed surface in the half-space z ≥ 0.

Taking this surface to be a half-circle (or half-sphere in 3D) of radius R joined with a segment of size 2R on the axis z = 0, we see that in the limit R → ∞

Φ(r 0 ) = C ∞ - κ D Sr dS(r s )Φ(r s )G N (r s |r 0 ), (24)
where we have used Eqs. ( 17),( 21),( 19) and ( 14) to simplify the integrals over the surfaces in Eq. ( 23). The above equation means that Φ(r) can be constructed for any position as soon as one knows its value on the reactive patch. Taking r 0 to be on the reactive patch yields an integral equation for Φ s , defined to be the value of Φ on the reactive patch. Since the above equation involves an unknown constant C ∞ it must be accompanied by a supplementary condition, which is provided by the relation (13). Let us finally write explicitly the integral equations for Φ for the 2D and the 3D cases. For d = 2, we obtain:

Φ s (x 0 ) = C ∞ + κ Dπ a -a dx Φ s (x) ln |x -x 0 |, (25) a -a dx Φ s (x) = 1/κ, (26)
In the case d = 3, we note that Φ s (r) depends only on the radial distance to the disk center, Φ s (r) = Φ s (r). The kernel of the integral equation can be simplified after a few algebraic manipulations detailed in Appendix B, leading to

Φ s (r 0 ) = C ∞ - κ D a 0 dr 2 r Φ s (r) π(r + r 0 ) K 2 √ rr 0 r + r 0 , (27)
2π a 0 dr r Φ s (r) = 1/κ, (28)
where K(k) is the complete elliptic integral of the first kind, defined as

K(k) = 1 0 dt[(1 -t 2 )(1 -k 2 t 2 )
] -1/2 with k the Elliptic modulus (to be distinguished from the parameter m = k 2 ). These integral equations admit no known analytical solution in general. In the next sections, we focus on their asymptotic study. From now on, without loss of generality, we set the units of length and time so that a = 1 and D/a 2 = 1. The remaining parameter κ then represents κa/D in full units.

III. THE LIMIT OF LARGE REACTIVITY

A. A scaling argument for the anomalous behavior of the mean reaction time for large reactivity.

Here we present a brief scaling argument that leads to the anomalous scaling (1). In the case of perfect reactions κ = ∞, it is clear from Eqs. ( 12),( 14),( 15) that Φ can be seen as the electrostatic potential generated by a charged conducting disk that is embedded in an insulating surface, with the prescription that Φ = 0 on this disk. In fact, using the image methods, it is easy to show that Φ is also the electrostatic potential in infinite space generated by an infinitely thin disk. It is well known [START_REF] Jackson | Classical electrodynamics[END_REF] that, for this electrostatic problem, the "electric field" -∇Φ diverges near the edge as 1/ρ 1/2 , where ρ is the distance from the disk extremity. Therefore, at a distance ρ a from the edges, Φ ∝ ρ 1/2 . In the case that κ is large but finite, since the natural length scale associated to finite reactivity is * = D/κ [START_REF] Grebenkov | Chemical Kinetics: Beyond the Text[END_REF], we may therefore assume that the mean reaction time is comparable to the mean first passage time when the starting position is located at a distance * from the reactive patch. In this condition, with ρ = * we obtain the anomalous scaling Φ ∝ 1/ √ κ, as announced in Eq. [START_REF] Redner | A guide to First-Passage Processes[END_REF]. In what follows, we show how to rigorously derive this scaling law, with the prefactor.

B. 2D case

Leading order

Consider now the limit κ → ∞ in the case d = 2. Since this situation corresponds to a first passage problem, we know that the asymptotic value of C ∞ does not depend on κ. The fact that Φ s vanishes in the large κ limit leads us to postulate in line with Eqs. ( 25)-( 26

) that Φ s (x 0 ) ∼ κ→∞ 1 κ Φ 1 (x 0 ), C ∞ ∼ κ→∞ C 1 , (29)
where C 1 and Φ 1 do not depend on κ. Inserting these ansatz into the integral equation ( 25) and the normalisation condition (26), we obtain

1 -1 dx Φ 1 (x) ln |x -x 0 | = -πC 1 , (30)
1

-1 dx Φ 1 (x) = 1. ( 31
)
The integral equation ( 30) for Φ 1 (x) is known as Carleman's equation and its analytical solution is known explicitly [START_REF] Polyanin | Handbook of integral equations[END_REF]. Using also the normalisation condition, we obtain the final expression for C 1 and Φ 1 :

C 1 = ln 2 π , Φ 1 (x) = 1 π √ 1 -x 2 . ( 32
)
2. Boundary layer near the extremities of the reactive patch

We now note that Φ 1 (x) is formally infinite at x = ±1, i.e. near the boundary of the reactive patch. This means that our expansion ( 29) is not valid near these points, suggesting a behavior similar to those obtained for boundary layer problems. Since the reaction length [START_REF] Grebenkov | Chemical Kinetics: Beyond the Text[END_REF] is 1/κ in our units, we expect processes happening at such scales. Therefore, we assume the behavior

Φ s (x) = κ α ψ((1 -|x|)κ), (1 -|x|) 1, ( 33
)
with ψ a scaling function. The exponent α will be set such that the behavior of Φ in the boundary layer matches with that far from the boundary layer. Namely, the compatibility of the above ansatz with Eq. ( 32) imposes the choice

α = -1/2, ψ(X) ∼ X→∞ 1 π √ 2X . (34)
Hence, the structure of the solution in the limit κ → ∞ is

Φ s (x) = κ -1 Φ 1 (x) + ... (1 -|x|) 1/κ κ -1/2 ψ((1 -|x|)κ) + ... (1 -|x|) 1 ( 35
)
and the condition (34) ensures that these two expressions give the same result in their common validity regime

κ -1 1 -|x| 1.
A key point here is that the mean return time, starting from the boundary of the reactive region scales as 1/κ 1/2 and is thus infinitely larger than the mean return time starting from the center of the target, which scales as 1/κ. This suggests to set

C ∞ = C 1 + C 1/2 /
√ κ + ... in the limit of large reactivity (even though the constant C 1/2 will turn out to vanish).

Let us now find the set of equations satisfied by ψ. First, we consider the normalisation condition (26), which we write under the form

1 -1 dx Φ s (x) - Φ 1 (x) κ = 0. ( 36
)
Here we remark that the integrand in the above integral is maximal near x = ±1. Let us define an intermediate length scale such that

1/κ 1. ( 37
)
We will keep this notation in the whole paper. We write Eq. ( 36) by separating the integral into two regions: when (1 -|x|) > then we approximate Φ s (x) by the first line of Eq. ( 35) and if (1 -|x|) < we use the expressions on the second line of Eq. ( 35) for Φ s and we approximate by its behavior near x = ±1, which reads

Φ 1 1/[π 2(1 -|x|)]. This leads to 0 = -1+ -1 dx ψ((1 -|x|)κ) √ κ - 1 κπ 2(1 -|x|) + 1 1- dx ψ((1 -|x|)κ) √ κ - 1 κπ 2(1 -|x|) .
Setting X = (1 -|x|)κ, and using κ 1, we obtain

∞ 0 dX ψ(X) - 1 π √ 2X = 0. ( 38
)
In order to find the equation satisfied by ψ it is useful to write the difference between the general equation ( 25) and the equation satisfied by Φ 1 , to find

Φ s (x 0 ) = C ∞ -C 1 + κ π 1 -1 dx Φ s (x) - Φ 1 (x) κ ln |x -x 0 |. ( 39
)
This leads to

Φ s (x 0 ) = C ∞ -C 1 + κ π -1+ -1 dx ψ((1 -|x|)κ) √ κ - Φ 1 (x) κ ln |x -x 0 | + κ π 1 1- dx ψ((1 -|x|)κ) √ κ - Φ 1 (x) κ ln |x -x 0 |. ( 40
) For x 0 = 1 -X 0 /κ, if we set X = (1 -|x|)κ, this yields ψ(X 0 ) = C 1/2 + κ 0 dX π ψ(X) - 1 π √ 2X ln |X -X 0 | κ + 1 π κ 0 dX ψ(X) - 1 π √ 2X ln 2 - X + X 0 κ . (41)
Now, in the limit κ → ∞, noting that κ → ∞ [by definition of , see Eq. ( 37)] and using the previously found condition (38) we obtain

ψ(X 0 ) = C 1/2 + ∞ 0 dX ψ(X) - 1 π √ 2X ln |X -X 0 | π . (42)
This is a Wiener-Hopf integral equation for the unknown function ψ(X). We solve it in the next section.

3. Solution of the Wiener-Hopf equation [START_REF] Deegan | [END_REF] We solve the Wiener-Hopf integral equation with Carleman's method, as described in Ref. [START_REF] Polyanin | Handbook of integral equations[END_REF]. We note that a similar equation has appeared in viscous flow theory [START_REF] Brown | [END_REF][47][48] but the differences between our equations and the equation studied in these references justify the fact to solve it here in detail. First, let us introduce the following notations. We denote f + (X) all functions (depending on the real variable X) that vanish for all X < 0, and f -(X) all functions that vanish for X > 0. For any function f (X) one can write

f (X) = f + (X) + f -(X), with f + (X) = f (X)θ(X) and f -(X) = f (X)θ(-X)
, where θ is the Heaviside step function. We introduce the complex Fourier transform and its inverse:

f (z) = ∞ -∞ dXf (X)e -izX , (43)
f (X) = 1 2π ∞ -∞ du f (u)e +iuX (44)
where z represents a complex number and u a real number. We denote f+ (z) the Fourier transform of the function f + (X), and f-(z) the Fourier transform of f -(z).

Typically, Fourier transforms of the form f+ (z) are defined in the lower complex half-plane Im(z) ≤ 0, Fourier transforms of the form f-(z) are defined in the upper complex half-plane Im(z) ≥ 0 (as long as f ± (x) does not diverge exponentially at x → ±∞). Now, we can define ψ + (X) ≡ ψ(X)θ(X), and we introduce

K(X) = 1 π ln |X|, χ + (X) = 1 π √ 2X θ(X). ( 45
)
The integral equation ( 42) can be generalized for negative X 0 by writing

ψ + (X 0 ) = ∞ 0 dX [ψ + (X) -χ + (X)] K(X -X 0 ) + y -(X 0 ), ( 46
)
where the only remarkable property of y -(X 0 ) is that it vanishes for positive X 0 . Note that we have assumed that C 1/2 = 0, this will be justified at the end of the calculation by the fact that the obtained solution satisfies the normalization condition (38) for this value of C 1/2 .

Taking the Fourier transform of the above equation, we obtain

ψ+ (u) = [ ψ+ (u) -χ+ (u)] K(u) + ŷ-(u). ( 47
)
Calculating the Fourier transforms leads to

ψ+ (u) 1 + 1 |u| = - 1 -i sign(u) 2 √ π|u| 3/2 + ŷ-(u), (48)
where sign(u) = θ(u) -θ(-u) is the sign function. This equation can be considerably simplified by introducing an auxilliary function S -(X) defined by

S -(X) = θ(-X) 2 |X| π , Ŝ-(u) = 1 -i sign(u) 2 √ π |u| 3/2 , ( 49
)
so that the Wiener-Hopf equation can be written as

ψ+ (u) 1 + 1 |u| = -Ŝ-(u) + ŷ-(u) = f-(u), ( 50
)
where f-(u) is the Fourier transform of another unknown function f -(X), whose only remarkable property is to vanish for positive X. This kind of equations is known as a homogeneous Wiener-Hopf equation, the method to solve it consists in obtaining a factorization of the form ψ+ (u)ĝ + (u) = ĝ-(u). To this end, we write the Wiener-Hopf equation under the form

ψ+ (u)e Ŵ (u) = f-(u), (51) with Ŵ
(u) = ln[1 + 1/|u|], (52)
and we seek a factorization Ŵ (u) = Ŵ+ (u) + Ŵ-(u). This can be done by calculating its inverse Fourier transform:

W (X) = cos(X)[2 Si(|X|) -π] + π 2π |X| - Ci(|X|) sin(X) πX (53)
where Ci and Si are the integral cosine and integral sine functions

Ci(X) = - ∞ X dt cos(t) t , Si(X) = X 0 dt sin(t) t . (54)
A factorization may thus be obtained by setting x). Now, we write the equation (51) as

W (x) = W + (x) + W -(x), i.e W + (x) = W (x)θ(x) and W -(x) = W (x)θ(-
ψ+ (u)e Ŵ+(u) = f-(u)e -Ŵ-(u) . (55)
We are now in the favorable case: the terms on left-handside are analytic functions in the upper complex plane, those on the right are analytic in the lower complex plane except for one pole at z = 0, and these terms are equal on the real axis. According to the theorem of analytic continuation, combined with the Cauchy theorem, we conclude that both terms are equal to a constant plus a 1/z term on the whole complex plane [START_REF] Polyanin | Handbook of integral equations[END_REF]. We thus have

ψ+ (u)e Ŵ+(u) = a 0 + a 1 u , (56)
where the constants a 0 , a 1 will be found by requiring that ψ(X) is a solution to our problem. Since Ŵ+ (z) is defined on the lower complex plane, we may consider the above equations on the lower imaginary axis u = -is, in which case the above equality can be written in terms of Laplace transforms, with the usual notation f+

(s) = ∞ 0 dte -st f + (t) = f+ (-is): ψ+ (s) = a 0 + i a 1 s e -W + (s) . (57)
The Laplace transform W + (s) can be identified by calculating its derivative, i.e. the Laplace transform of -xW (x), and then by integrating over s; this leads to

W + (s) = 1 4 ln 1 + s 2 s 2 + m(s), (58) with m
(s) = - s 0 dw ln w π(1 + w 2 ) . (59)
We know that the behavior of ψ for large arguments is given by the matching condition Eq. ( 34), which translates to the small-s behavior:

ψ+ (s) ∼ s→0 1 √ 2πs . ( 60
)
Inserting Eq. ( 58) into (57) and taking the small-s limit, we see that the above behavior is obtained for ia 1 = 1/ √ 2π. Next, the value of a 0 is found by requiring that ψ(0) is finite, so that ψ(s) vanishes in the limit s → ∞; this leads to a 0 = 0. Hence, the final expression for the function ψ is given in the Laplace domain by

ψ+ (s) = 1 √ 2πs(1 + s 2 ) 1/4 e -m(s) . (61)
Finally, we must check that the normalization condition (38) holds, this is readily done by noting that ψ(s) -1/ √ 2πs = O( √ s ln s) vanishes for small s. This justifies our hypothesis C 1/2 = 0. Unfortunately, the Laplace inversion cannot be performed so that we know only ψ(X) in closed form, however we can easily derive the asymptotic behavior of ψ(X) for small and large arguments. The study of the asymptotic behavior of ψ(s) for large s leads readily to the initial value of ψ(X):

ψ(0) = lim s→∞ s ψ(s) = 1 √ 2π , (62)
this justifies the previously announced result [START_REF] Redner | A guide to First-Passage Processes[END_REF]. The large X behavior can be computed by noting that the Laplace transform of Xψ(X) is d ψ/ds and by expanding this one for small s, with the result:

ψ(X) X→∞ 1 √ 2Xπ + ln(4X) + γ e -1 π 2 (2X) 3/2 + ... (63)
with γ e the Euler-Mascheroni constant. A formula for ψ(X) can be obtained by considering the inverse Laplace transform of ψ(s) with the Mellin's inverse formula, by using a contour that follows the negative real axis (above and below), such Laplace inversion is obtained in Appendix 4 of Ref. [48]:

ψ(X) = 1 π ∞ 0 dp e -pX+m(p) √ 2πp(1 + p 2 ) 3/4 . ( 64
)
In summary, here we have obtained an analytic expression in Laplace space for the scaling function ψ(X) which characterizes the behavior of the mean first passage time near the extremities of the reactive patch in two dimensions. The validity of our approach is checked on Figure 3 by comparing with exact numerical results obtained from the general form of the solution.

FIG. 3. Behavior of the mean first reaction time near the extremities of the reactive patch. Symbols: exact general solution obtained numerically in Appendix A in 2D (upper symbols) and 3D (lower symbols). We also represent the values of ψ and ψ 3d obtained from Eqs. ( 64) and (78).

Next-to-leading order expansion

Up to now, the constant C ∞ , which characterizes the behavior of the mean reaction time when the initial position is far from the target, has been obtained at leading order only in the large reactivity limit, with the same result as in the case of a perfectly reactive patch. Here we show how to obtain the first non-trivial correction for C ∞ for large reactivity, with the result that C ∞ (κ) is not analytic in κ. First, we note that when 1/κ 1 -|x| 1, Eqs. ( 63) and (35)

indicate that Φ s (x) 1 κπ 2(1 -|x|) + ln[4(1 -|x|)κ] + γ e -1 κ 2 π 2 [2(1 -|x|)] 3/2 + ... (65)
This suggests that, outside the boundary layer, the nextto-leading order behavior of Φ s reads:

Φ s (x) = Φ 1 κ + Φ * 2 ln κ + Φ 2 κ 2 + ..., (66)
because this expression can be matched with (65) by imposing that

Φ * 2 (x) 1 π 2 [2(1 -|x|)] 3/2 , (x → ±1), ( 67
) Φ 2 (x) ln[4(1 -|x|)] + γ e -1 π 2 [2(1 -|x|)] 3/2 , (x → ±1). ( 68
)
These expansions also lead us to assume that, for large reactivity the constant C ∞ behaves as

C ∞ = C 1 + C 2 + C * 2 ln κ κ + ... (κ → ∞) (69)
The equation for Φ 2 and Φ * 2 can be identified as follows. We consider again the intermediate length scale satisfying (37), and we start from the integral equation ( 39), which we write as

1- -1+ dx Φ * 2 (x) ln κ + Φ 2 (x) κ ln |x -x 0 | = -π C 2 + C * 2 ln κ κ + πφ 1 (x 0 ) κ + B(x 0 ) + B(-x 0 ), ( 70
)
where B contains all the terms which appear due to the fact that the integral over Φ 2 in the above equation is evaluated over a truncated interval ] -

1 + ; 1 + [ instead of ] -1; 1[, so that B = κ 0 dX √ κ 1 π √ 2X -ψ(X) ln 1 -x 0 - X κ . ( 71
)
To proceed further, we consider (70) as an integral equation for Φ 2 + ln κΦ * 2 over the truncated interval ] -1 + ; 1 + [. Its solution is analytically known and we identify the constants C 2 and C * 2 by requiring that the normalisation condition is satisfied at this order of κ. This procedure requires to evaluate B in the limit → 0 without assuming that 1 -|x 0 |, and it turns out to be relatively tedious. The calculation is described in Appendix C 1, and leads to the explicit results

C * 2 = 1 π 2 , C 2 = 1 + γ e + ln 8 π 2 . ( 72
)
These values of C 2 and C * 2 are in excellent agreement with the exact solution for Φ s obtained numerically, as shown in Fig. 4(a).

C. 3D case

We now adapt the approach to the 3D case. It turns out that the solution admits the same scaling behaviors than in 2D:

Φ s (r) = Φ1(r) κ + Φ * 2 (r) ln κ+Φ2(r) κ 2 (1 -r) 1/κ, 1 √ κ ψ 3d ((1 -r)κ) (1 -r) 1, ( 73
)
where the first line is the expansion of Φ s (r) in powers of κ at fixed r, and the second line the expansion of Φ in powers of κ at fixed X = (1 -r) √ κ. At leading order, the integral equation for Φ 1 reads

0 = C 1 - 1 0 dr 2r π(r + r 0 ) K 2 √ rr 0 r + r 0 Φ 1 (r), ( 74
) 2π 1 0 dr r Φ 1 (r) = 1. ( 75
)
The solution of the above integral equation (where r Φ 1 (r) is considered to be the unknown function) is known [START_REF] Polyanin | Handbook of integral equations[END_REF] and this leads to the solution

Φ 1 (r) = 1 2π √ 1 -r 2 , C 1 = 1 4 . (76)
We thus note that

Φ 1 (r → 1) ∼ 1 2π 2(1 -r) . (77)
In the boundary layer near r = 1, we set r = 1 -X/κ, r 0 = 1 -X 0 /κ, and Φ s (r) = 1/ √ κψ(X). With these scalings we can expand the integral equation (27) with the result that ψ 3d satisfies exactly the same equation than in 2D, the only difference is that it has to match with ψ 3d (X) ∼ 1/(2π √ 2X) for large X [due to Eq. ( 77)] and there is thus a factor of 2 that arises when we compare to the 2D case:

ψ 3d (X) = 1 2 ψ(X). ( 78
)
This relation is checked on Fig. 3. Let us now identify the next order terms in 3D. Inserting the ansatz (73) into the integral equation ( 27) and expanding at second order, we obtain

1- 0 dr r r + r 0 K 2 √ rr 0 r + r 0 [Φ 2 (r) + Φ * 2 (r) ln κ] = - πΦ 1 2 + π 2 [C * 2 ln κ + C 2 ] + B ( 79
)
where the term B compensates the fact that the above integrals are evaluated over the truncated interval [0; 1 -[, so that:

B(r 0 , ) = - √ κ 1 1- dr r r + r 0 K 2 √ rr 0 r + r 0 × ψ 3d ((1 -r)κ) - 1 2π 2(1 -r)κ . (80)
As in the 2D situation, we consider (79) as an integral equation for which the solution is analytically known; and we then chose C 2 and C * 2 so that the normalisation condition for Φ s holds at all orders of κ. The final result is

C * 2 = 1 4π , C 2 = γ e + 1 + ln(2) 4π , (81)
and it agrees perfectly with numerical solutions, as shown in Fig. 4(b).

IV. THE LIMIT OF SMALL REACTIVITY

Let us now consider the limit κ → 0. At leading order, the mean reaction time is homogeneous. We seek a solution under the form

Φ s (x) = 1 κ n≥0 f n (x)κ n , C ∞ = 1 κ n≥0 c n κ n . (82)
At leading order, we obtain

f 0 = c 0 = 1/|S r |. ( 83
)
where S r is the length (in 2D) or the area (in 3D) of the reactive patch. Furthermore, next-orders can be found iteratively by using

f n (r) = c n - 1 D Sr dS(r )f n-1 (r )G N (r|r ) (84)
with the condition for n ≥ 1:

Sr dSf n = 0. ( 85
)
For d = 2, the explicit computations can be done for the first orders, and we find

c 0 = 1/2; c 1 = 3 -ln 4 2π , c 2 = 2 9 - 7 3π 2 . (86)
In the 3D situation, the leading order is simply

f 0 = c 0 = 1/π; (87)
and the recurrence relation is

f n (x) = c n - 2 π 1 0 dyf n-1 (y) y x + y K 2 √ xy x + y . (88)
Unfortunately, it seems very difficult to calculate these integrals, and even at first order the coefficient c 1 can be calculated only numerically:

c 1 = 4 π 2 1 0 dx 1 0 dy xy x + y K 2 √ xy x + y 0.27. ( 89
) V. COMPARISON WITH THE CONSTANT FLUX APPROXIMATION
The constant flux approximation (CFA) [20] has been used in many recent studies [37,39,40] on imperfect reactivity in confinement, and here we consider how this approximation compares to the exact results in our formalism. First, we need to adapt this approximation to our situation of large volume limit. In the CFA, one replaces the Robin condition ( 14) by inhomogeneous Neumann conditions:

D∂ z Φ = 0 (z = 0, r > a) -Q (z = 0, r < a) , (90)
where the flux Q is assumed to be constant on the reactive patch and will be determined self-consistently with a closure relation. A natural choice of closure relation is to impose that Robin condition is satisfied on average, hence

Q = κ Sr dS Φ s , (91)
but we also have the normalization condition (13), so that

Q = 1. (92)
Now, inserting (90) into ( 23) leads directly to a solution for Φ within CFA:

Φ(r 0 ) = C ∞ + Q D Sr dS(r s )G N (r s |r 0 ). (93)
Integrating over S and using (92) and (91), we obtain the CFA value of C ∞ :

C cfa ∞ = 1 κS r 1 - κ D Sr dS(r) Sr dS(r 0 )G N (r|r 0 ) (94)
Comparing with the results of Sec. IV, we see that in the CFA, C ∞ is exactly the same as the next-to-leading order expansion of C ∞ in the limit of low reactivity, i.e.,

C cfa ∞ = c 0 κ + c 1 (95)
It may be therefore surprising that CFA works for C ∞ even for rather large values of the reactivity (Fig. 5), but this comes from the fact that the value of c 1 turns out to be extremely close to the exact value of C ∞ (κ = ∞) (the difference is of the order of a few percents). This might be the reason why the CFA approach can be implemented to yield accurate results in other contexts. However, the value of the mean first passage time near the extremity of the reactive patch is not well captured by this approximation, since it is obvious in Eq. ( 93) that it does not scale as 1/ √ κ for large κ, contrary to what we have found. FIG. 5. Values of C∞ in 2D (a) and 3D (b), as found from the exact numerical solution compared to the exact large and small reactivity asymptotics. Note that the constant flux approximation (CFA) is exactly equivalent to the first order expansion in the limit of low reactivity. Here we use the units so that D = 1, a = 1.

VI. CONCLUSION

In this paper we have considered the imperfect narrow escape problem for diffusive particles in confinement. We have established a general formalism which provides the mean reaction time in the large volume limit for any value of the reactivity parameter. We have obtained explicit results in d = 2 and d = 3 in the respective limits of low and large reactivity parameter. Our most surprising result is the scaling of the mean reaction time when the initial position is at the extremity of the imperfect patch; this mean return time scales as κ -1/2 and is thus much larger than the naively expected scaling 1/κ. Interestingly, we have shown that this anomalous scaling is closely related to the divergence of the electric field near corners and edges of conducting objects [START_REF] Jackson | Classical electrodynamics[END_REF], which is also responsible for the existence of coffee rings [START_REF] Deegan | [END_REF] or the crispiness of the extremities of cooked potatoes [49]. We have explicitly identified the prefactor of this scaling law by solving a Wiener-Hopf equation. We have also identified a non-analytic behavior for the capacitances of the imperfect patches as a function of the reactivity. We note that we have restricted ourselves to the case of circular patches, but we believe that for the more general patches with a smooth boundary the asymptotic scaling laws should remain unchanged. Finally, we have made a link between the results obtained within the Constant Flux Approximation (CFA) and the low reactivity limit. It turns out that the CFA gives very accurate predictions of the mean reaction time when the initial position is far from the target, but fails to predict the correct behavior of the mean return times; when the initial position is on the reactive patch. In the future, one could adapt our formalism to multiple targets, for example to generalize the classical calculation [12] of the absorption time by a sphere covered by reactive patches to imperfect patches. Our results provide a general framework to quantify the mean reaction times for the imperfect narrow escape problem. with (r, θ, ϕ) the usual spherical coordinates. It is useful to calculate the scale factors h i = |∂r/∂i|, with i = {α, β, ϕ},

h α = a α 2 + β 2 1 + α 2 1/2 , h β = a α 2 + β 2 1 -β 2 1/2 , h ϕ = a (1 + α 2 )(1 -β 2 ) 1/2 . ( A16
)
For axisymmetric functions, the Laplacian reads in this orthogonal coordinates

∇ 2 Φ = 1 h α h β h ϕ ∂ ∂α h β h ϕ h α ∂Φ ∂α + ∂ ∂β h α h ϕ h β ∂Φ ∂β , , (A17)
so that Φ satisfies the equation

∂ ∂α (1 + α 2 ) ∂Φ ∂α + ∂ ∂β (1 -β ) ∂Φ ∂β = 0. (A18)
We impose Neumann conditions for β = 0 and β = 1, at which ∂ β Φ = 0. With these conditions, the general solution can be found by the method of separation of variables, which leads to

Φ(α, β) = Φ ∞ + ∞ q=0 a q g q (α) P 2q (β), (A19)
where P 2q are even Legendre polynomials (satisfying both Neumann conditions at β = 0 and β = 1), and

g q (α) = 1 i Q 2q (iα) - π 2 P 2q (iα), ( A20
)
where i 2 = -1 and Q 2q are Legendre functions of the second kind. Let us give here additional details on the function g q . To see that g q is real it is useful to write Q 2q as [START_REF] Abramowitz | Handbook of Mathematical Functions With formulas, Graphs, and Mathematica Tables[END_REF] Q

2q (x) = P 2q (x) 2 ln 1 + x 1 -x -W 2q-1 (x), ( A21
)
where W is the polynomial

W 2q-1 (x) = q m=1 4q -(1 + 4(m -1)) (2m -1)(2q -m + 1) P 2q-(2m-1) (x). (A22)
For purely imaginary arguments x = iα, we have

Q 2q (iα) = iP 2q (iα) arctan(x) -W 2q-1 (iα), (A23)
and we thus see that

g q (α) = P 2q (iα) arctan(α) - W 2q-1 (iα) i - π 2 P 2q (iα). (A24)
Using the parity of P and W , it becomes clear that g q is real. Furthermore it can be checked that it decreases to zero at infinity (and g 0 ∼ 1/α for large α).

Now, the equation satisfied by the coefficients a q is identified by using the Robin condition. In these coordinates, the partially absorbing disk corresponds to α = 0, and the Robin conditions can be deduced from

∂ n Φ = -h -1 α ∂ α Φ α=0 so that the boundary conditions read (D∂ α Φ -aβκΦ) α=0 = 0. (A25)
Inserting the general solution (A19) into the above boundary condition, multiplying by P 2k (β) and integrating, we obtain the linear system: ∞ q=0 a q g q (0)

1 0 dβP 2q (β)P 2k (β) = κΦ ∞ 1 0 dββP 2k (β) + κ ∞ q=0 a q g q (0) 1 0 dββP 2q (β)P 2k (β), (A26)
for all positive integers k. Finally, we calculate the surface element at α = 0, dS α = h ϕ h β dϕdβ so that the normalization condition reads

2πκa 2 1 0 Φ(0, β)βdβ = 1, (A27)
which leads to the equation

κa 2 πΦ ∞ + 2πκa 2 ∞ q=0 a q g q (0) 1 0 dββP 2q (β) = 1. (A28)
A numerical solution for Φ ∞ can thus be found by solving the linear system (A26) for the coefficients a q (completed by the above normalisation condition). Note also that C ∞ = Φ ∞ .

Appendix B: Identification of the integral equation ( 27) in 3D

Here we briefly show how to obtain the integral equation (27). Using Eq. ( 20) for d = 3, we see that Eq. ( 24) writes

Φ(r 0 ) = C ∞ - κ D a 0 drK(r, r 0 )Φ(r), (B1) with K(r, r 0 ) = 1 2π 2π 0 r dθ r 2 + r 2 0 -2rr 0 cos θ . ( B2
)
The quantity K(r, r 0 ) can be simplified as follows:

K(r, r 0 ) = 1 π π 0 r dθ r 2 + r 2 0 -2rr 0 cos θ = 2r π π/2 0 du r 2 + r 2 0 -2rr 0 cos(2u) = 2r π π/2 0 du r 2 + r 2 0 -2rr 0 [2 cos 2 u -1] = 2r π(r + r 0 ) π/2 0 du 1 -4rr0 (r+r0) 2 cos 2 u = 2r π(r + r 0 ) K 2 √ rr 0 r + r 0 , (B3)
where u = θ/2 and in the last line we have recognized the definition of the elliptic function K. Inserting the above result into Eq. (B1) finally leads to Eq. ( 27).

B - 1 √ κ κ 0 dX ψ(X) - 1 π √ 2X ln(1 -x 0 ) - κ 0 dX √ κ ψ(X) - 1 π √ 2X ln 1 -x 0 -X/κ 1 -x 0 . (C1)
Here we have only assumed that one can use the leading order of the scaling form for the solution Φ for arguments lower than . Now, to evaluate the first line we use the normalisation condition (38), and to evaluate the terms on the second line, we change variable X = uκ :

B 1 √ κ ∞ κ dX ψ(X) - 1 π √ 2X ln(1 -x 0 ) - √ κ 1 0 du ψ(κu ) - 1 π √ 2κu ln 1 -x 0 -u 1 -x 0 . (C2)
Using Eq. ( 63), we can write B under the form

B f * 2 (x 0 ) ln κ + f 2 (x 0 ) κ , (C3) with f * 2 (x 0 ) = ln(1 -x 0 ) √ √ 2π 2 - 1 0 du π 2 (2u ) 3 2 ln 1 -x 0 -u 1 -x 0 , (C4) f 2 (x 0 ) = ln(4 ) + γ e + 1 √ 2 π 2 ln(1 -x 0 ) - 1 0 du ln(4u ) + γ e -1 π 2 (2u ) 3/2 ln 1 -x 0 -u 1 -x 0 . (C5)
Let us precise a few properties of f * 2 (similar properties hold for f 2 ). In the limit → 0 at fixed x 0 , we see that

f * 2 (x 0 ) ∼ →0 ln(1 -x 0 ) √ √ 2π 2 (C6)
Near the extremity of the patch, we set x 0 = 1 -v to determine the behavior of f * 2 . In the limit → 0 at fixed v = (1 -x 0 )/ , we obtain We consider this equation as an integral equation over the interval [-1+ ; 1-], for which the solution is analytically known [START_REF] Polyanin | Handbook of integral equations[END_REF]: When → 0, the behavior of I * 2 ( ) is obtained by inserting the small limit of f * 2 (x 0 ) at fixed x 0 given by Eq. (C6) into Eq. (C9) and inserting the result into the above equation, leading to

f * 2 (x 0 = 1 -v ) ∼ →0 ln(v ) √ √ 2π 2 - 1 0 du π 2 (2u )
Φ * 2 (x) = 1 π 2 √ b 2 -x 2 b -b dt √ b 2 -t 2 ∂ t F * 2 (t) t -x + 1 ln[b/2] b -b dt F * 2 (t) √ b 2 -t 2 , (
I * 2 ( ) →0 (-1) π ln 2 1 -1 dt ln 1 -t 2 √ 2π 2 √ √ 1 -t 2 = √ 2 π 2 √ , (C12)
this result is consistent with the fact that the behavior of Φ * 2 near x = ±1 is given by (67). Now, the fact that the normalisation condition (26) ) .

(C14)

The contributions of f * 2 (t) in this integral are negligible except for x 0 at the vicinity of 1. Thus, we set t = 1 -v to evaluate the above integral, and the integral can be evaluated by integrating v over [0, ∞] (except for the term C * 2 which is evaluated without this change of variable). Using (C7) yields

∆I * 2 ( ) (-1) π ln 2 1 -1 dt -πC * 2 √ 1 -t 2 + 2 π 2 √ 2 ∞ 0 dv ln(2v ) θ(v -1) 2(v -1) - 1 √ 2v - ∞ 1 dv 1 0 du 2 π 2 (2u) 3/2 2(v -1) ln v -u v . (C15)
All the integrals appearing in the above equations can be evaluated. Imposing ∆I * 2 ( ) = 0 then leads to

C * 2 = 1/π 2 . ( C16
)
To identify C 2 we proceed in the same way. The integral equation is

1- -1+ dx Φ 2 (x) ln |x -x 0 | = F 2 (x), (C17)
F 2 (x) = -πC 2 + πφ 1 + f 2 (x 0 ) + f 2 (-x 0 ), (C18) so that

I 2 ( ) = b -b dx Φ 2 (x) = 1 π ln(b/2) b -b dt F 2 (t) √ b 2 -t 2 . (C19)
As before we need to evaluate the behavior of F 2 (x, ) for small

F 2 (x) ∼ →0 ln(4 ) + γ + 1 √ 2 π 2 ln(1 -x 2 0 ) = F 0 2 (x) √ . (C20)
At leading order for small we obtain

I 2 ( ) ∼ →0 1 π ln(1/2) 1 -1 dt F 0 2 (t) √ √ 1 -t 2 = I 0 2 √ . (C21)
The next-to-leading order is

I 2 - I 0 2 √ (-1) π ln 2 1 -1 dt F 2 (t)θ(b -|t|) √ b 2 -t 2 - F 0 2 (t) (1 -t 2 ) (C22)
Again, we evaluate it by setting t = 1 -v and taking the small limit of the obtained integrand at fixed v, leading to To evaluate the term involving Φ 1 we introduce a variable ε so that ε 1 and we calculate b

-b dx Φ 1 (x) √ b 2 -x 2 = b -b dx π √ 1 -x 2 √ b 2 -x 2 1-ε -1+ε dx π(1 -x 2 ) + ε/ 1 2 dv 2π v(v -1 ln(8/ ) π , (C24)
where the last equality follows from the evaluation of the integrals with ε 1. Finally, evaluating all terms in Eq. (C23) and requiring that the I 2 ( ) -I 0 2 ( ) → 0 for small , we obtain

C 2 = 1 + γ e + ln 8 π 2 , ( C25
)
as announced in the main text.

2. Second-order calculation in the limit of large reactivity in 3D

We evaluate the term B in Eq. (80) by writing )

B(r 0 , ) √ κ K 2 √ r0 1+r0 (1 + r 0 ) ∞ κ dX ψ(X) - 1 2π √ 2X - 1 0 du √ κ(1 + r 0 ) ψ(u κ) - 1 2π √ 2u κ × K 2 (1 -u )r 0 1 -u + r 0 -K 2 √ r 0 1 + r 0 . ( C26

FIG. 1 .

 1 FIG. 1. (a) Illustration of the imperfect narrow escape problem. A partially reactive patch (thick red line) is embedded in the boundary of a confining domain. A random walker, starting from the initial position (red sphere) diffuses in the domain and is eventually adsorbed on the patch. (b) Zoom on the portion of space delimited by the dashed blue lines around the reactive patch.

  FIG. 2.(a) Geometry of the confining domains (called A and B) that are considered for stochastic simulations. In 2D, these domains are defined in polar coordinates by r(θ) = Rf (θ) with f = 1.6(1 + 0.5 cos 2 θ) for domain A and f = 1.6(1+0.1 sin θ+0.3 sin 3θ) for domain B. Domains in 3D are obtained by considering revolution of 2D curves around the vertical dashed line. The reactive patch is indicated by a thick red line, and the initial position is taken at a distance r from the center of the patch along the black dashed line. In the figure, we have used R = 6a. (b),(c) Results of Brownian dynamics simulations for the mean reaction time in 2D/3D (parameters are indicated in the legend) compared to general theoretical expressions as obtained in Appendix A. In all simulations, we used a time step ∆t = 10 -4 a 2 /D. Boundary conditions are implemented as follows: if, at the end of a time step, the random walker falls outside the domain, then if it is "behind" a reflecting wall it is reflected with respect to this wall, and if it falls "behind" the reactive patch, it is absorbed with probability Pa = κ πdt/D (in which case the trajectory ends) and it is reflected with probability 1 -Pa, see Ref.[43].

FIG. 4 .

 4 FIG. 4. Comparison of the values of C∞ obtained numerically (symbols, see Appendix A for details), with the analytical predictions in Eqs. (72) and (81) (black lines), for the twodimensional (a) and three-dimensional (b) domains.

2 1

 2 Appendix C: Calculation details for the identification of C2 and C * . 2D situation Here we describe the details of calculations leading to the identification of the constants C 2 and C * 2 in the 2D situation. Let us first evaluate the term B in Eq. (71):

  O(ln κ/κ) in the integral equation (70) leads to1--1+ dx Φ * 2 (x) ln |x -x 0 | = F * 2 (x 0 ),(C8)F * 2 (x 0 ) = -πC * 2 + f * 2 (x 0 ) + f * 2 (-x 0 ). (C9)

  C10) where we have set b = 1 -. As a consequence, the integral of Φ * 2 over the truncated interval ] -b; b[ reads

I 2 (γ e - 1 π 2

 212 (2u) 3/2 2(v -1) ln v -u v (C23)

  holds at all powers of κ

	implies that
			lim →0	I * 2 ( ) -	√ π 2 √ = 0. 2	(C13)
	We thus evaluate
	∆I * 2 ( ) = I * 2 ( ) -	√ π 2 √ 2
	(-1) π ln 2	1 -1	dt	F * 2 (t)θ(b -|t|) √ b 2 -t 2	-	ln 1 -t 2 π 2 2 (1 -t 2

ACKNOWLEDGMENTS

Computer time for this study was provided by the facilities MCIA (Mesocentre de Calcul Intensif Aquitain) of the Université de Bordeaux and of the Université de Pau et des Pays de l'Adour. T. G. acknowledges the support of the grant ComplexEncounters, ANR-21-CE30-0020.

Appendix A: Exact general form of the solution Φ(r)

Imperfect narrow escape problem in 2D

Here we describe a way to obtain the general solution of the problem formed by Eqs. ( 12)- (14). It consists in writing the equations in a set of orthogonal coordinates and using the standard method of separation of variables. We first describe this approach in the 2D case, for which we use the elliptic coordinates µ, ν defined as x = a ch(µ) cos(ν), y = a sh(µ) sin(ν).

(A1)

We calculate the scale factors h i = |∂ i r| with i ∈ {µ, ν}:

The Laplace equation satisfied by Φ and the reflecting boundary conditions outside the reactive patch are written in this coordinate system as

The general solution for these equations can be written by using the method of separation of variables, which leads to

φ n e -2nµ cos(2nν). (A4)

Furthermore, the normalization condition (13) can also be written D n Φ = 1 for any surface surrounding the target. Far from the target, this means that ∂ r Φ = 1/(πrD). Noting that µ ln(2r/a) for large r, we thus find

We also note that the quantity C ∞ is given, in this mode decomposition, by

Finally, the Robin condition at the target surface reads

Using this condition and the form of the general solution (A4), we find that the coefficients φ n are solution of the infinite linear system

which is satisfied for all positive integers m, with

In practice, this linear system (A9) can be solved numerically by taking into account only a finite number of modes N , and checking that the obtained quantities do not depend on N for large N . Note also that C ∞ can be directly calculated by using Eq. (A6).

3D case

This approach can be adapted to the 3D case, for which we use orthogonal coordinates defined as

where ϕ is the azimuthal angle. Note that α > 0 and β ∈ [0; 1] are related to the standard oblate spheroidal coordinates (µ, ν, ϕ) by α = sinh(µ) and β = sin(ν). Inversion formulas read

Note that here, for conciseness we will treat ln κ as being of order 1 in powers of κ, the result will be exactly the same as in the case where one separates the ln κ terms and the O(1) terms.

In the small limit at fixed r 0 we obtain

whereas if we set r 0 = 1 -v, in the limit → 0 at fixed v we obtain

where we have used K(1 -y)

2 ln(8/y) for small y. Let us write the integral equation (79) under the form

Using the analytically known solution [START_REF] Polyanin | Handbook of integral equations[END_REF] of the integral equation (C30), we obtain

When → 0 we obtain at leading order

and this integral diverges for → 0, as it should due to the known behavior for Φ2 (r) when r approaches 1. At next-to-leading order, we evaluate the terms involving B by setting s = 1 -v and take the small limit at fixed v, so that we can use Eq. (C29):

To evaluate the term containing Φ 1 , defined in Eq. ( 76) we can use again a trick where we use a variable ε with ε 1: Finally, all the integrals in (C34) can be evaluated, requiring that it vanishes for small leads to

which is exactly Eq. (81).