Kinetic modeling of protein PEGylation

Author's Accepted Manuscript

Kinetic Modeling of protein PEGylation David Pfister, Eric Bourgeaux, Massimo Morbidelli

www.elsevier.com/locate/ces

PII: DOI: Reference:

S0009-2509(15)00516-3 http://dx.doi.org/10.1016/j.ces.2015.07.031 CES12503

To appear in:

Chemical Engineering Science

Received date: 5 February 2015 Revised date: 4 June 2015 Accepted date: 27 July 2015 Cite this article as: David Pfister, Eric Bourgeaux, Massimo Morbidelli, Kinetic Modeling of protein PEGylation, Chemical Engineering Science, http://dx.doi.org/ 10.1016/j.ces.2015.07.031 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting galley proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Kinetic Modeling of Protein PEGylation David Pfister†, Eric Bourgeaux‡, Massimo Morbidelli†,* †

Institute for Chemical and Bioengineering, ETH, Vladimir-Prelog-weg 1, 8093 Zurich,

Switzerland ‡

DSP Laboratory, Sanofi, 9 quai Jules Guesde, 94400 Vitry-sur-Seine cedex, France

*Corresponding author: phone: +41 44 632 30 34; fax: +41 44 632 10 82. E-mail address: [email protected] (M. Morbidelli).

Abstract PEGylation is currently the technique of choice to extend the in-vivo circulation half-life of therapeutic proteins. Nevertheless, the PEGylation reaction is not yet fully described, particularly concerning the development of a proper kinetic model. In this work, we propose a detailed model which describes the rate of protein amine PEGylation with functionalized PEG esters by accounting for the specific reactivity of the various residues and for the shielding effect generated by the already conjugated PEG chains. The proposed kinetic rate constants are derived so as to explicitly include diffusion limitations, heterogeneous reactivities and the steric hindrance caused by the polymer layer surrounding the PEGylated proteins. With illustrative purposes, this new model has been applied to the simulation of lysozyme PEGylation at the lysine residues with PEG esters of various molecular weights and with various PEG to protein molar ratios. Keywords: PEGylation; Kinetic Modeling; Positional Isomers; Isoforms Distribution

1

1.

Introduction

Numerous techniques have been developed over the past decades to enhance the in-vivo stability of therapeutic proteins (Birch and Onakunle, 2005; Saltzman, 2001). In particular, PEGylation, known as the covalent binding of PEG chains to proteins, has become a technique of choice to extend their in-vivo circulation half-life (Abuchowski et al., 1977a; Abuchowski et al., 1977b; Pfister and Morbidelli, 2014; Veronese et al., 2009; Zalipsky and Harris, 1997). Despite its great popularity, PEGylation suffers from several drawbacks. Indeed, the additional costs caused by the expensive PEG functionalization, the conjugation reaction and the supplementary purification steps dramatically impact the cost of the final product. The implementation of kinetic models can bring into prospective competitive advantages when it comes to efficiently develop and design a process while minimizing the associated experimental work and development cost. Indeed, the development of a kinetic model can provide a deeper understanding of the experimental observations and can help for the optimization of the operating conditions. There is often much more information to extract from experimental results by means of a kinetic model than those obtained by simple inspection of the raw data. Moreover, in the frame of the Quality by Design (QbD) initiative (ICH Q8 (ICH, 2005a, 2007), Q9 (ICH, 2005b) and recently Q11 (ICH, 2012)), biopharmaceutical industries must adhere to very strict norms in terms of quality and safety to ensure reproducibility and robustness of their production. Therefore, QbD principles turn model-based approaches to powerful tools to identify the critical parameters affecting the process. Laura et al. (Laura et al., 2007) were the first to propose a simple kinetic model to simulate Growth Hormone Antagonist (GHA) PEGylation. Their kinetic model is easy to implement and represents an effective tool to simulate the formation of PEGylated proteins at various degrees of PEGylation (also named PEGamers). This model has already been used with some success 2

in the literature (Molek, 2008; Moosmann et al., 2011; Wang et al., 2014). The main limitation of the approach proposed by Laura et al. is that the employed kinetic model does not allow the distinction between the various positional isomers that compose PEGamers of a given degree of PEGylation (doP). Since positional isomers (also named isoforms) are known to exhibit different reactivities (Maiser et al., 2012), this lumped strategy might lead to strong discrepancies with the experimental data. The model developed in this work intends to give a better insight into the PEGylation kinetics by accounting for all the positional isomers. By avoiding the lumping procedure, the full distribution of isoforms is obtained which may become in the future a requirement for defining the quality of the final product. Moreover, contrary to the existing kinetic model, the kinetic rate constants determined in this work explicitly account for the two important consequences of the PEGylation reaction that affect the protein reactivity: an increase of the hydrodynamic radius which leads to a reduced diffusion rate and the shielding effect caused by the conjugated PEG polymers. For the sake of generality, the kinetic model has been developed for any kind of reactant exhibiting in general

N conjugation sites. To do so, each isoform is described by a binary

coding indexation. This new mathematical treatment, which is here used to describe the mass balance of each isoform during the PEGylation reaction, can also be used for modeling any kind of reaction involving molecules with multiple reacting sites, such as protein-inhibitor binding, selective olefin functionalization, multiple protecting group addition, partial hydrogenation of oil, etc. The predictions of the detailed PEGylation model are compared to experimental data obtained with various PEG to protein ratios and PEG molecular weights.

3

2.

Model Development

2.1.

Kinetic scheme

The kinetic model developed in this work provides a detailed description of the PEGylation process without resorting to any lumping strategy. The model accounts explicitly for two different types of reaction: -

a bimolecular stepwise reaction between the functionalized polymer and a conjugation site of the considered protein;

-

in parallel, the reaction of hydrolysis of the active moiety of the functionalized PEG.

The hydrolysis reaction of the PEG moiety is modelled as a pseudo first-order reaction, with effective kinetic rate constant kd . For the sake of generality, a generic protein with N different residues suitable for conjugation is considered. The protein is therefore discretized with respect to its number of PEGylation sites. In the following example each index corresponds to a site ( N in total) whose value is 0 for a free site or 1 for the already PEGylated site. The stepwise stoichiometric conjugation of this protein with the functionalized PEG can be written as follows:

k

1

k

2

Protein 0,0,…,0,0 + PEG 0→ Protein1,0,…,0,0 Protein 0,0,…,0,0 + PEG 0 → Protein 0,1,…,0,0 ... k

i

Protein 0,0,…,0,0 + PEG 0→ Protein 0,0,.,1,.,0,0

where PEG and Protein are the concentration of functionalized PEG and protein, respectively. Later on, the concentration of the native protein will be noted P0 while Pj will be the PEGamer concentration of doP j .

4

The conjugated proteins are described both in terms of doP (number of conjugated PEG) and in terms of isoforms distribution for a given doP. The number of possible isoforms can be computed from the doP j and the total number of conjugation sites N as follows: N ! (N − j ) ! j ! . The total number of species is 2N which leads to complex and computationally demanding models when N is increasing.

In order to generalize the kinetic model to any kind of protein with N PEGylation sites, the involved species were differentiated by means of a binary coding indexation. For instance, the binary coding of a protein having 3 PEGylation sites is the following: -

000 corresponds to the native protein;  

-

100 , 010 and 001 are the 3 different mono-PEGylated positional isomers;      

-

110 , 101 and 011 are the 3 different di-PEGylated positional isomers;      

-

and finally, 111 corresponds to the fully PEGylated protein.

The corresponding distribution of the different species is shown in Figure 1a. The binary vectors mentioned above are then compiled in a ( 2N , N ) matrix B whose lines correspond to the binary coding of integers going from 0 (for the native protein) to 2N − 1 (for the fully PEGylated protein).  0 1 0 1 0 1 01    B =  0 0 1 1 0 0 11   0 0 0 0 1 1 11    T

The advantage of this method is to be adaptable to any kind of protein since it is generalizable to N sites. However, since it considers all the possible positional isomers, it requires also a

tremendous number of kinetic rate constants. Considering N conjugation sites would in fact require N × 2N + 1 kinetic rate constants for the ‘full model’ instead of the only N + 1 needed in the existing model. Such detailed description of the protein reactivity is only possible when the number of sites is limited, as for cysteine PEGylation for instance. In the case where a large 5

number of sites is involved, even though the mathematical description of the ‘full model’ would be possible, analytical methods would usually fail to identify all the species required for evaluating the large number of model parameters. This difficulty can be resolved by introducing the following fundamental assumption: -

the kinetic rate constant of a given site is only affected by the doP and not by the position of the already conjugated PEGs on the protein.

This means that the conjugation rate constant on a specific site for all proteins Pj containing j PEG polymers, is the same for all the positional isomers having the same value of j . This assumption leads to the ‘site-specific reactivity’ model since the heterogeneity in the reactivity of the PEGylation sites is assumed to be independent of the location of the already conjugated polymers. Under this assumption, the number of kinetic rate constants can be reduced down to N × N + 1 . This is not to be confused with the existing model where, in addition, the different sites

are assumed to react one after the other, reducing the number of independent rate constants to N + 1 as indicated above. The model proposed by Laura et al. could be named ‘infinite-

selectivity’ model since each PEGamer is assumed to be composed of a single positional isomer, thus reducing the distribution of isoforms to a single element. There is therefore a perfect correspondence between PEGamer and isoform in their model.

6

Figure 1: a) Distribution of positional isomers of a protein with 3 PEGylation sites (
: free site,  PEGylated site). b) Kinetic rate constants associated with different models (red: full model, blue: site-specific reactivity model and black: iso-reactivity model). In order to illustrate the model proposed in this work, let us consider the reaction of conjugation of a PEG molecule on the third position of a mono-PEGylated protein described by the vector 010 , i.e. the already conjugated PEG is located on the second position of a protein having   3 three position sites. The corresponding kinetic rate constant is denoted k010 , where the exponent

stands for the targeted PEGylation site, while the subscript describes the PEGylation state of each site. By assuming that only the number of conjugated PEG polymers (but not their location) affects the kinetic rate constants of a given reactive site, the notation introduced above can be 7

simplified: within the above mentioned assumption, the subscript can be replaced by the doP. 3 simplifies to k13 . Note Therefore, continuing the previous example, the kinetic rate constant k010

that this is different from the iso-reactivity assumption discussed later for which all the sites are assumed to have the same reactivity which in this case would be denoted k1 (for which there is no need for a subscript any longer). This notation is illustrated in Figure 1b.

2.2.

Kinetic model

Using the formalism previously introduced, it is possible to generalize the mass balance for all the different species, from the native to the fully conjugated protein, as follows:

dci dt

(

)

(

)

= (δi ⋅ B) ⋅ kn −1
ci* PEG − δi ⋅ ( 1 − B) ⋅ kn ci PEG i

i

(1)

where ci , i = 0, 2N − 1 , is the molar concentration of species i whose description in terms of doP and conjugated sites is given by the line i + 1 of the matrix B (the vector containing all the ci is noted c ), δi is the Kronecker row vector whose elements equal 0 except the i + 1 coordinate which equals 1, and
is the Hadamard product (element-by-element) of two vectors. n is the vector of degrees of PEGylation obtained by summing the elements of each row of B . ci* is the vector of N species leading to the formation of ci , 1

is the ( 2N , N ) matrix whose elements

equal 1, and k is a column vector of N kinetic rate constants whose dependency on the doP ni will be described in the next section. ci* is computed from the binary coding of i but with free sites on the positions 1 to N , computed

as follows:

8

 c   i −1   c   i −2   c   i −4  * l −1 ci =   if i − 2 ≥ 1, l ∈ [1, N ] (0 otherwise) c  i −8     ...    c N −1   i −2 

(2)

Continuing with the same example considered above, the formation of the di-PEGylated protein conjugated at the positions 2 and 3 can only derive from the mono-PEGylated protein conjugated at the position 2 (with no PEG on the site 3) and from the mono-PEGylated protein conjugated at the position 3 (with no PEG on the site 2). In other words, the binary vector 011 (decimal value: i =6) derives from 001 (decimal value 4) and 010 (decimal value 2). According to

equation

(2)

c*6 = c5 c4 c2  ,

( δ6 B ) . ( k
c*6 ) = k11 × 0 + k12 × c4

T

and

δ6 B =  011  ,

so

it

comes

+ k13 × c2 which only contains the two species involved in the

formation of the considered isoform 011 . Two additional equations must be computed to account for the consumption of PEG both by conjugation with proteins and hydrolysis. The high reactivity of the functionalized PEG goes in fact along with the hydrolysis of its active moiety which cannot be neglected. Moosmann et al. (Moosmann et al., 2011) showed that accurate simulations are not possible without considering this deactivation reaction. Esters hydrolysis mechanisms have been thoroughly investigated over the years and appear to be pH-dependent, affected by both basic and acidic conditions (Day and Ingold, 1941). Generally esters hydrolysis kinetics are assume to be pseudo-first order reactions while the kinetic rate constant accounts for all deactivation sources, mostly pH. Under this simplifying assumption the functionalized PEG hydrolysis rate is described by the following first order kinetic reaction: dPEGd dt

= kd PEG

9

(3)

where kd is the pseudo-first order kinetic rate constant of the functionalized PEG deactivation and PEGd is the hydrolysed PEG concentration. Finally, for a closed batch reactor, the total mass of polymer is conserved. Consequently the rate of consumption of the functionalized polymer is given by: dPEG dc = −kd PEG − nT . dt dt

(4)

Thus, the total number of ordinary differential equations (ODE) to be solved simultaneously is 2N + 2 .

2.3.

Kinetic rate constants

2.3.1. Diffusion-limited kinetic model When dealing with macromolecules, as polymers and proteins, the diffusion of the reactants starts playing a role in the dynamics of reaction. For not so rapid chemical reactions compared to the approach of the reacting species, the characteristic time of diffusion becomes comparable to the characteristic time of reaction. Therefore, the kinetic rate constants result from the contribution of two effects in series: the rate of approach and the rate of the subsequent chemical reaction. Moreover, the rate of the chemical reaction is hindered by the so-called ‘shielding effect’ which is a consequence of the presence of a PEG layer surrounding PEGylated proteins. Laura et al. (Laura et al., 2007) reported for example that very high degrees of PEGylation are never observed, even at high PEG/protein ratios, due to steric hindrances and repulsive interactions between the free functionalized PEG in solution and the PEG layer surrounding the conjugated protein. The rate of approach has been derived based on the pioneering work of von Smoluchowski (von Smoluchowski, 1906, 1917). Since then, many efforts have been put forth to develop mechanistic models to describe diffusion-limited chemical reactions (Calef and Deutch, 1983; 10

Collins and Kimball, 1949; Hammes, 1978; Noyes, 1961, 1969). The use of a diffusion-limited kinetic model was motivated by the experimental observation that the protein PEGylation rate is decreasing when increasing the polymer molecular weight and the doP. In order to give rise to the chemical reaction step, the two reactants (protein and PEG) have to diffuse toward each other while interacting in a potential field U( r ) , deriving from the repulsive/attractive forces between PEG and protein. The potential strength depends on the centre-to-centre distance r between the two molecules. The probability to find a PEG molecule in the neighbourhood of a PEGylated protein at a distance r is described by the distribution function p (r, t ) which is governed by the diffusion equation. At infinite distance p (r, t ) tends to the bulk value given by the product of the concentrations of the two species, Pj × PEG (Hammes, 1978; Schulten, 2000). Otherwise, the variations of p (r, t ) in time and space are described by the following equation:

∂p = DPEG + DP j ∂t

(

)

 ∂2p 1 ∂Fp    −  ∂r 2 kBT ∂r 

(5)

where p (r, t ) represents the average concentration of reactants pairs separated by a distance r and F is the sum of the forces acting on the two molecules deriving from the potential energy ∂U ∂r = −F . DPEG and DP

j

are the diffusion coefficients of the PEG and conjugated protein,

respectively. In order to derive the diffusion-limited kinetic rate constants, the diffusion equation (5) has been solved under stationary conditions. Considering that the chemical reaction only occurs at the contact point r0 = RPEG + RP0 at the surface of the native protein, where RPEG and RP

0

are the hydrodynamic radii of the PEG and native protein respectively, the boundary

condition at r0 is described by the flux equation:

11

Jtot (r0 ) 4πr02 = hp(r0 )

(6)

for a given constant h describing the effectiveness of the chemical reaction. h , in m/s, can be seen as a mass transfer coefficient; the larger h , the faster the chemical reaction. Jtot is the flux of pairs (protein, PEG) coming into contact and therefore represents the total reaction rate. Recalling that ∂p ∂r − pF kBT = e

(

−U kBT

U kBT

∂ e

p

) ∂r , the total flux estimated at the distance r

is given as follows:

J tot (r ) = 4πr

2

(D

PEG

)

+ DP e j

−U kBT

(

U kBT

∂ e

p

)

(7)

∂r

The conservation of the total flux induces that Jtot is the same for all r and must hold in particular at the contact distance r0 . Moreover the potential is assumed to be nil at infinite distance, so that the integration of p (r ) between r0 and infinity leads to:

p ( r0 ) =

(D (D

PEG

+ DP

PEG

+ DP

j

j

)+r

2

0

) P PEGe

he

j

−U 0 / kBT

−U 0 / kBT ∞ eU (r )/ kBT

∫r

0

r2

(8) dr

The integral of the potential is estimated by the Fuchs’ approximation (Fuchs, 1934): ∞

r0 ∫ r0

U k BT

e

r

2

U m k BT

dr ≈ λe

(9)

where λ is a constant deriving from the integration of the l.h.s of equation (9) and Um is the extreme value of the potential energy profile. According to Yang et al. (Yang et al., 2011) the PEG/protein interaction is attractive, resulting into a negative interaction potential. For PEG and insulin, they estimated this potential to be about -60 kJ/mol. Finally, conventionally defining the rate of bimolecular conjugation as k j Pj PEG , it follows that:

12

kj =

(

)(

)

c RPEG + RPc DPEG + DP 4π −U m / kBT j 0 e c λ 1 + κ DPEG + DP RPEG + RPc

(

j

)(

0

(10)

)

where κ equals e(U 0 −Um )/kBT λh . The term at the denominator, κ ( DPEG + DPj

) (R

c PEG

+ RPc

0

)

is

the inverse Thiele modulus of the process which represents the ratio of the characteristic time of diffusion over the characteristic time of reaction. Compared to equation (6), the radii for the calculation of the contact distance were replaced by the collision radii Rc = 1.65R according to (Sandkühler, 2004). In order to be able to apply the model given by equation (10) for the PEGylation rate constants, we need to compute the interaction potential, U , and the diffusion coefficients of the two reacting species , DPEG and DP , and finally to quantify the shielding effect of the conjugated PEG chains j

on the PEGylated proteins reactivity.

2.3.2. Site Specific reaction It is known that the reactivity of each protein residue depends on its nucleophilicity, surface accessibility and pKa (Lee and Richards, 1971; Suckau et al., 1992). The reactivity of amine groups of lysozyme was investigated toward acetylation (Suckau et al., 1992), PEGylation (Lee and Park, 2003) and estrone glucuronide acylation (Smales et al., 1999) and in all cases the most reactive residues were Lys97 and Lys33, which corresponds to the highest accessibility values. The reactivity of each amine is therefore driven by the surrounding residues which affect the interaction potential of the lysine residues of interest. With a view to developing a model as general as possible, site specific reactivity has been implemented in equation (10) by setting specific values of the interaction potential for each site. In other words: U m ≡ U mi . Moreover, the parameter U 0 , which is the value of the interaction potential at contact, is also specific for each site. In the following, it will be assumed that the variations of U 0,i among the 13

isoforms are such that the difference U 0,i −U m,i could be considered as a constant for all the residues. Consequently, an average value for κ is used in equation (10).

2.3.3. Diffusion coefficients In the frame of the Stokes-Einstein theory (Einstein, 1905; von Smoluchowski, 1906) the diffusivity of spherical particles is estimated as follows:

D=

kBT

(11)

6πηR

The hydrodynamic radii (in Å) are calculated according to the core-shell model introduced by Fee and Van Alstine (Fee and Van Alstine, 2004) in which both PEG and protein are considered as spheres. Under this assumption, the hydrodynamic radius is directly linked to the molecular weight of each species by a power law. RP = 0.82 ( M wProt 0

)

R P E G = 0.19 ( M wP E G

1/ 3

(12)

)

(13)

0.56

The hydrodynamic radius of the PEGylated protein, RPj , is calculated as follows:

RP = j

A 2 2 1 + RPEG ( j ) + RPEG ( j ) 6 3A 3

  3 3 (j ) + 12 81RP6 + 12RP3 RPEG (j ) with A =  108RP0 + 8RPEG 0 0  

(

1 1 3 2

)  

(14)

and where RPEG ( j ) is the radius of a

hypothetic polymer whose molecular weight is equal to the sum of the already conjugated polymer chains ( = jM wPEG ).

14

2.3.4. Shielding effect Once attached to the protein, the polymer evolves around the protein and forms an obstacle that covers partially the surface of the PEGylated protein (Fee and Van Alstine, 2004; Meng et al., 2012). Moreover, molecular dynamic simulations showed that conjugation with a PEG chain results in a shielding effect which reduces the solvent accessible surface area (SASA) of the protein. Therefore, the shielding effect lowers the reactivity of the PEGylated protein by masking free sites and prevents them from being conjugated. According to MD simulations with low molecular weight polymers, the SASA has a logarithmic decay with the PEG molecular weight (Yang et al., 2011). Consequently, motivated by the work of Yang et al., it is assumed that the fraction of accessible surface area can be expressed as follows:

 SASA   = −αjMwPEG ln   SASA0 

(15)

where the parameter α is proportional to the average surface of the protein covered by a monomer unit. A more detailed description of the ‘shielding effect’ requires extensive MD simulations which is not the purpose of this work. Therefore, the empirical equation (15) will be considered later on as it provides a simple quantification of this effect with respect to the PEG molecular weight. Thus, in conclusion, substituting diffusivity coefficients by equations (11) in equation (10) and multiplying the kinetic rate constant by the fraction of accessible surface area leads to the kinetic rate constant of PEGylation on the site i of a protein already conjugated with j PEG molecules:  2k T i k j i =  B e −U m / kBT  3λη

(

)(

)

 c RPEG + RPc 1 RPEG + 1 RP   PEG j 0   e −α jM w     1 + κ ' 1 R c c RPEG + RP 1 PEG + 1 RPj  0

(

)(

)

    2

(16)

which can simply be written as k ji = k0i ϕ j with k 0i = { ... }1 and ϕ j = {... }2 . In this frame, the total number of kinetic parameters is given by: 15

-

N ‘intrinsic’ kinetic rate constants k 0i for the lysine residues;

-

one hydrolysis rate constant kd ;

-

and two parameters α and κ ' = κ × kBT 6πη .

It is worth noting that in the case of chemically-controlled reactions, implying large κ ' values,

(

PEG

−α jM equation (16) simplifies to k j i = 6.6πr02e w

)

λκ e

i −Um /kBT

, since the effect of diffusion is

negligible. One easily recognise an Arrhenius-like kinetic rate constant. These kinetic rate constants are consistent with concentrations expressed in number per volume. For molar based concentrations k0i should be multiplied by the Avogadro constant.

3.

Materials and methods

3.1.

Chemicals

All the PEGylation experiments considered in this work were performed in 25 mM phosphate buffer at pH 7.0 by mixing a 6 g/L lysozyme (HEWL, 98% pure from Sigma) solution with mPEGSPA (JenKem), with nominal average molecular weight of 5 kDa (PEG5), 10 kDa (PEG10) and 20 kDa (PEG20). The buffer solutions for the chromatographic experiments were prepared with analytical grade chemicals. The water was deionized and further purified with a Millipore Simpack 2 purification pack (Millipore). All the buffers were filtered on Durapore membrane filters, type HVLP, 0.45 µm (Millipore). The chemicals for SDS-PAGE and the protein standards were purchased from Life Technologies.

16

3.2.

Experiments

3.2.1. Reaction procedure 0.5, 1, 1.5, 3 and 6-fold molar excess of mPEG-SPA were dissolved in 6 g/L lysozyme solution at pH 7.0. The reaction took place at 25ºC in a stirred beaker during 3h. Samples were taken at different times and quenched by rapid addition of 0.1 M HCl to drop the pH below 3.0 (Molek, 2008).

3.2.2. PEG characterization 3.2.2.1. Molecular weight MALDI-TOF analysis was used to determine the exact molecular weight of the mPEG-SPA. The measurements were performed with MALDI-TOF Ultraflex II (Bruker Daltonics). The exact molecular weights were 5.4, 10.6 and 20.2 kDa, respectively. 3.2.2.2. Hydrolysis rate 2 mg of mPEG-SPA were dissolved in 1.0 mL of 25 mM phosphate buffer adjusted to pH 7.0. The rate of hydrolysis was monitored by UV spectroscopy at 260 nm by observing the formation of free N-hydroxysuccinimide (NHS) from the cleavage of the ester bond between the propionic acid and the NHS. The experiments were performed at room temperature. The kinetic rate constant

of

hydrolysis,

kd

was

determined

from

the

linearized

absorption

plot:

− ln ( ( A∞ − At ) ( A∞ − A0 ) ) vs time. The slope directly corresponds to the hydrolysis rate

constant (Harris and Kozlowski, 1997). The kinetic analysis has been performed with 5, 10 and 20 kDa mPEG-SPA. Results are shown in Figure 2, where it can be seen that the PEG molecular weight has no influence on the rate of hydrolysis.

17

Figure 2: Kinetic analysis of mPEG-SPA hydrolysis monitored by UV spectroscopy at 260 nm. Normalized absorbance of the free NHS group released during hydrolysis of 5 kDa mPEG-SPA at room temperature in 25 mM phosphate buffer at pH 7.00. In the inset are presented the linearized data for 5, 10 and 20 kDa mPEG-SPA. kd , obtained by linear regression is constant within the different investigated molecular weights and equal to 9.7 10-3 min-1.

3.2.3. Separation of the PEGamers The characterization of the PEGylated lysozyme was performed on an Agilent Series 1100 unit equipped with a diode-array detector. The peaks were detected at 280 nm. 3.2.3.1. Size Exclusion Chromatography. Crude solution from the reaction was analysed by SEC on an analytical TricornTM Superdex 200 10/300 GL (GE Healthcare). As mobile phase a 25 mM sodium phosphate buffer, containing 100 mM Na2SO4 at pH 7.0 was used. The elution flow rate was set to 0.5 mL/min. 3.2.3.2. Ion-Exchange Chromatography The cation-exchange resin Gigacap S-650M resin (Tosoh Bioscience) was chosen for the semipreparative purification of mono-PEGylated lysozyme. The stationary phase was packed in a Tricorn column (GE Healthcare) according to the manufacturer protocol. The resulting bed length was 38 mm. The mobile phase was a 20 mM citrate buffer at pH 3.0. The buffers were prepared with citric acid (Sigma-Aldrich) for buffer A while buffer B was also containing 1 M sodium 18

chloride (Millipore). The buffer solutions were adjusted to pH 3.0 by addition of sodium hydroxide as needed and filtered. The volumetric flowrate was set to 1 mL/min. The column was preequilibrated with 15% B. A piecewise gradient strategy was applied: 26% B over 18 minutes, 60% over 5 additional minutes and 100% B over the last 5 minutes. 100% B was maintained during 4 minutes and the column was reequilibrated at 15% B during 6 minutes. The pure monoPEGylated lysozyme was collected between 20 and 25 minutes. The fractions were then pooled together. If needed, pooled fractions were concentrated with Vivaspin Turbo 15 (Sartorius) during 1 hour at 4000 rcf.

3.2.4. Characterization of the positional isomers 3.2.4.1. pH-gradient elution chromatography A Dionex Propac SCX-10 analytical column (4x250 mm) preceded by the corresponding guard column (4x50 mm) from Thermo Fischer Scientific was used for the separation of the positional isomers. The separation was performed by pH-gradient elution chromatography according to Maiser et al. (Maiser et al., 2012). 20 mM CABS (Sigma) buffer solutions adjusted to 10.50 and 11.50 for buffer A and B respectively were used. The pH was monitored with a FTPH micro flow through pH meter from LASAR Research previously calibrated with high precision standard solutions. The elution was carried out with a linear gradient from 0% to 100% B over 15 minutes, then 5 minutes at 100% B and the column was reequilibrated during 10 minutes at 0% B. The flow rate was set to 1 mL/min. 3.2.4.2. SDS-PAGE For the SDS-PAGE experiments NuPAGE® Bis-Tris mini gels were used and the analysis was carried out in XCell SureLock® Mini-Cell. The gels were prepared according to the manufacturer instructions and run at constant voltage (200V) for 35 minutes. The upper and lower chambers were filled with 1xNuPAGE MES running buffer. Concentrated samples from the crude and purified mono-PEGylated protein were diluted with deionized water. All the samples were mixed with NuPAGE LDS sample buffer (1/4 in volume) and heated during 10 minutes at 70°C. The gel 19

was then treated according to the SilverXpress® silver staining protocol. The protein standards was Mark12™ Unstained Standard from Life Technologies, with insulin A and B chains (2.5 and 3.5 kDa), Aprotinin (6.0), Lysozyme (14.4), Trypsin inhibitor (21.5), Carbonic anhydrase (31.0), Lactate dehydrogenase (36.5), Glutamic dehydrogenase (55.4), BSA (66.3), Phosphorylase b (97.4), β-galactosidase (116.3) and Myosin (200).

3.3.

Numerical methods

The model developed in this work was implemented in Matlab®. The fitting parameters ( k0i , α and κ ' ) were obtained by minimizing the residual sum of squares between experimental and simulated distribution of isoforms and PEGamers kinetics using the embedded function fminsearch.

4.

Results and discussion

4.1.

Kinetic analysis

As mentioned above, the kinetics of PEGamers formation have been monitored by SEC. An example of the obtained chromatograms is shown in Figure 3 for the case of a 6 g/L lysozyme solution PEGylated with 1.5-fold molar excess of PEG5. The samples were taken at 1, 3, 5, 10, 20, 30, 45, 60, 90, 120 and 180 minutes.

20

Figure 3: Kinetic analysis of the PEGylation reaction with 6 g/L lysozyme solution and a PEG to protein molar ratio of 1.5. The evolution of the concentrations was followed by SEC and monitored by UV absorbance at 280 nm. The peaks have been assigned by comparing the retention times to standards. The last eluting compound is very small and corresponds to the NHS fragment released by conjugation or hydrolysis of the functionalized PEG. Then, the different species elutes according to their molecular weight starting from the more retained lysozyme to the early eluting tetra-PEGylated protein. No doP higher than 4 has been observed. With four PEG chains conjugated at the surface of the protein the shielding effect is probably so pronounced that no additional PEG can reach the surface of the protein. From the data shown in Figure 3, it is seen that lysozyme is consumed while the doP progressively increases in time. When the concentration of lysozyme becomes sufficiently small, the concentration of mono-PEGylated protein also starts decreasing, and so on, until no functionalized PEG is left in solution. In order to obtain the concentration profiles in time of the different species, the peaks were numerically deconvoluted. The analysis was performed in Origin® by assuming a gaussian shape for all the peaks. The results of such a deconvolution process are presented in Figure 4 for three chromatograms obtained at certain specific operating conditions for PEG5, PEG10 and PEG20, respectively. The peaks deconvolution method is a powerful technique to analyse poorly 21

resolved peaks inevitably observed when increasing the doP and the PEG molecular weight. Finally, as the PEG is invisible in UV at 280 nm, its concentration could only be estimated by integrating equation (4) and closing the mass balance.

Figure 4: Example of peaks deconvolution. The experimental data (
) correspond to 6 g/L lysozyme PEGylated with 1.5-fold molar excess of PEG5 and 3-fold molar excess of PEG10 and PEG20. The individual gaussian curves are depicted as dashed lines while the continuous line correspond to the sum of all the signals.

4.2.

Positional isomers distribution

The key parameters for the detailed kinetic model are the ‘intrinsic’ kinetic rate constants of each site, k0i . Their proper estimation obviously requires the measurement of the distribution of the positional isomers of mono-PEGylated lysozyme as a function of time. For this, the protocol developed by Maiser et al. (Maiser et al., 2012) was applied. Typical chromatograms are shown in Figure 5. Accordingly, the mono-PEGylated lysozyme was firstly fractionated from the native protein and multi-PEGylated proteins. Based on previous work by Moosmann et al. (Moosmann et al., 2010), a separation procedure by IEC was developed. As suggested by Moosmann (Moosmann et al., 2010) and Maiser et al. (Maiser et al., 2012), the cation exchange resin 22

Toyopearl S-650M from Tosoh Bioscience was selected. Among the different investigated pH conditions (data not shown), pH 3.0 gave the best resolution between mono- and di-PEGylated lysozyme. The gradient strategy was optimized to get the best resolution in the shortest time while keeping the mono-PEGylated protein at the highest possible concentration. Fractions were collected, their purity assessed by SEC, and - if pure enough - were pooled together. The separation was repeated for the samples taken at 1, 5, 20, 45 and 180 minutes. A typical chromatogram corresponding to the PEGylation of lysozyme with 1.5-fold molar excess of PEG5 is shown in Figure 5a. It can be seen that the highly PEGylated species were not bound under these conditions and eluted first.

Figure 5: (a) Chromatogram of the separation of lysozyme PEGamers by cation-exchange chromatography on Toyopearl S-650M. The doP was confirmed by SEC and comparison with standards proteins. (b) Resolution of the mono-PEGylated lysozyme positional isomers by pH 23

gradient elution chromatography on PROPAC SCX-10. Numbers 1 to 6 corresponds to the different positional isomers. The pooled fractions were further analysed by pH gradient elution chromatography on a strong cation-exchange resin. By using a linear ascending pH gradient crossing the isoelectric point (pI) of the different mono-PEGylated molecules (between 10.5 and 11), Maiser was able to achieve superior resolution for the separation of the isoforms compared to classical salt gradient (Maiser et al., 2012). As expected, the pI of the mono-PEGylated proteins is slightly lower than the pI of the native lysozyme (11.2) due to the consumption of amine residues and to the shielding effect. Therefore, a linear gradient between pH 10.5 and 11.5 was applied. CABS, 4(Cyclohexylamino)-1-butanesulfonic acid, was the salt of choice for this high pH buffer solution, and the selected column was a Dionex PROPAC SCX 10 analytical column. The chromatogram of the separation of the isoforms from the above mentioned sample is presented in Figure 5b. The peaks were identified by SDS-PAGE. The results of the SDS-PAGE together with the corresponding chromatogram are shown in Figure 6. Five clear peaks were obtained, eluting after 10.6, 14.3, 15.4, 16.0 and 16.9 minutes, respectively. In addition, two partially resolved peaks at 11.2 and 14.7 minutes form shoulders after the two first well defined peaks. To confirm the nature of each peak the chromatogram was fractionated every 30 seconds. The fractions 1 to 7, as shown in Figure 6, were analysed by SDS-PAGE and compared with the purified monoPEGylated lysozyme, the crude solution from the reaction and to protein standards.

24

Figure 6: Chromatographic separation of the positional isomers of mono-PEGylated lysozyme resolved by pH gradient elution chromatography. The doP of the samples 1 to 7 was analysed by SDS PAGE. The main peaks were confirmed to correspond to mono-PEGylated lysozyme as well as the first barely resolved peak eluting after 11.2 minutes. On the contrary, traces of multi-PEGylated lysozyme were observed in fraction 4 suggesting that this peak corresponds to residual multiPEGylated lysozyme. Although lysozyme has 7 PEGylation sites (one N-terminal and six εamine of lysine residues), Maiser suggested that there is no possibility to distinguish the α from the ε-amine of Lys1. These six identified PEGylation sites are shown in Figure 7. In other words, each peak would correspond to one residue so that all the reactive sites are identified.

25

Figure 7: Lysozyme with its 6 PEGylation sites. ε-amine of Lys1 and N-terminal are lumped in a single site. The analyses were repeated with samples taken after 1, 5, 20, 45 and 180 minutes. As an example, the chromatograms corresponding to mono-PEGylated lysozyme isoforms obtained with 1.5-fold molar excess of PEG5 are shown in Figure 8. The chromatograms are plotted both with respect to elution time and reaction time. From these, the concentration of each isoform is computed and its evolution in time is also shown in Figure 8. It can be seen that the reactivities of the six sites are quite different and in particular two distinct groups can be observed: one composed by the three most reactive amine residues and the other composed of three less reactive amines. This is coherent with the findings of Lee et al. (Lee and Park, 2003), who identified Lys97, Lys33 and Lys116 as the most reactive residues towards amine PEGylation. The last 3 peaks were therefore attributed to the less reactive residues: Lys1, Lys13 and Lys96.

26

Figure 8: Kinetics of PEGylation of the mono-PEGylated lysozyme positional isomers. The analysed samples were taken at 1, 5, 20, 45 and 180 minutes, quenched in acidic solution, and separated by IEC on Gigacap S-650M. The mono-PEGylated thus obtained was analysed by pH gradient elution chromatography on PROPAC SCX 10.

4.3.

Kinetic Parameters Estimation and Model Simulations

As mentioned above, the developed kinetic model includes the hydrolysis kinetic rate constant, kd and the PEGylation rate constants, k ji which, using equation (16), involves in turn six ‘intrinsic’

rate constants k0i with i = 1, 6 and the two parameters α and κ ' . The first of the above parameters, kd was estimated by linearizing the experimental data presented in Figure 2 which led to kd = 9.7 10−3 min-1. It is worth noting that this value is comparable with the one found by Molek by direct fitting of PEGamers distribution data: kd = 1.0 10−2 min-1 (Molek, 2008). The ‘intrinsic’ kinetic rate constants, k0i have been estimated from the distribution of positional isomers shown in Figure 9. The shielding effect parameter, α as well as the parameter κ ' were fitted simultaneously on the distribution of PEGamers.

27

A comparison between fitted and experimental results is shown in Figure 9b. The experimental data correspond to the PEGylation of lysozyme with 1.5-fold molar excess of PEG5. The best fitted values of the ‘intrinsic’ rate constants, k0i with i = 1, 6 are 19.0, 1.9, 4.0, 3.7, 20.8 and 17.4 mL.mol-1.min-1, respectively. The so-estimated kinetic rate constants were used unchanged to predict the evolution of the distribution of the mono-PEGylated isoforms with only 0.5-fold molar excess of PEG5. Figure 9a shows the corresponding comparison between the model predictions and the experimental data.

Figure 9: Comparison between calculated and experimental mono-PEGylated lysozyme positional isomers distributions in the case of: (a) PEGylation with 0.5-fold molar excess PEG5, and (b) PEGylation with 1.5-fold molar excess PEG5. Solid curves are the model results and the points of the measured data of the six positional isomers: (
) [100000], () [010000], () [001000], () [000100], () [000010], and () [000001]. Kinetic parameters values obtained by fitting the PEGylation data with 1.5-fold molar excess PEG5 were used to predict the data for 0.5fold molar excess PEG5. The parameters α and κ ' were fitted to 4.42 10-5 mol.g-1 and 3.15 103 Å2, respectively. The values of the parameters used to simulate lysozyme PEGylation under nine different experimental conditions (three PEG/protein ratios and three PEG molecular weights) are summarized in Table 1. 28

Table 1: Model parameters Parameters Value -3

Fitting data

-1

kd

9.7 10 min

k 01

19.0 mL.mol-1.min-1

k 02

1.9 mL.mol-1.min-1

k 03

4.0 mL.mol-1.min-1

k 04

3.7 mL.mol-1.min-1

k 05

20.8 mL.mol-1.min-1

k 06

17.4 mL.mol-1.min-1

α

κ'

UV spectroscopy

4.4 10-5 mol.g-1 3.2 103 Å2

Positional isomers distribution

PEGamers distribution

In conclusion, the results of the simulations are compared with all considered experimental data in Figure 10. It appears that the model predictions are in fairly good agreement with these experimental data. In particular, it is seen that the model properly accounts for the effects of the PEG molecular weight (diffusivity and shielding effect) on the PEGylation kinetics up to a PEG molecular weight of 10 kDa. The data obtained with the largest PEG (20 kDa) show some discrepancy with the model prediction indicating we are certainly reaching the limit of the validity of equation (15). With this kinetic model in hands, the operational parameters (i.e. the PEG to protein molar ratio and the reaction time) can be optimized to produce PEGylated protein of a given doP at the maximum yield. Moreover, according to equation (16), the impact of the PEG molecular weight on the kinetic rate constants can be quantified a priori. The knowledge of the reaction kinetics is of primary importance to design and scale up efficient PEGylation processes. The detailed description of the protein reactive sites enables to use this model with any kind of PEGylation chemistry. For instance, functionalized PEG aldehyde are known to be more specific towards N-terminal PEGylation, but not entirely selective. The reactivity profile of the residues would certainly present a strong heterogeneity with a high reactivity value for the N-terminal residues and lower reactivities for the lysine residues. In the case of even more selective 29

chemistries as with cysteine PEGylation the number of sites N often resumes to one single site. The detailed description of the isoforms is consequently not necessary but the kinetic rate constant can still be described in the frame of the present model to account for the effect of the PEG molecular weight on the kinetics of PEGylation. To be complete, it should also be noted that some functionalized PEG are less sensible to hydrolysis and equation (4) can be modified so as to neglect this reaction.

Figure 10: Comparison between simulations (solid lines) and experimental data for different PEG concentrations and 3 different PEG molecular weights ((
) P0, () P1, () P2, () P3, () P4, and () PEG).

30

4.4.

The ‘iso-reactivity’ Model

The kinetic model developed in the previous section requires the determination of the reactivity of each residue. In the absence of proper analytical method, the distribution of isoforms required to estimate the specific kinetic rate constants cannot be obtained. Nevertheless, the kinetic model presented here can still be used at the expense of the detailed isoform distribution description. In the latter it is assumed that all active sites on the protein have the same reactivity, that is the ‘intrinsic’ rate constants k 0i = k for all i = 1, N . This is not to be confused with the existing model developed by Laura et al. (Laura et al., 2007) which assumes ‘infinite-selectivity’, since the residues are assumed to react one after the other. In order to simplify the detailed model presented in this work, the concentration of PEGamer is explicitly written as follows: 2N

Pl =

∑δ ( l − n j )c j

(17)

j =1

 0 if j ≠ l . Pl is the total concentration of conjugated protein with a doP l (sum of  1 if j = l 

with δ(l − j ) = 

the isoforms of doP l ). It is worth noting that equation (1) could be written in terms of pseudocomponent P0 to PN by considering the concentration dependent kinetic rate constants: N

k0 =

∑

k 0i i =1

 2N  c and k j = ∑k j  ∑δ ( j − nl )( 1 − Bli ) l P  l =1 j i =1 N

i

  .  

As a consequence of the ‘iso-reactivity’ assumption, the mole fraction of isoforms for a given doP are all equal. Under this simplifying assumption the concentration-dependent kinetic rate constant simplifies to:  N − 1  N     = Nk ( 1 − j N ) ϕj k j = N k ϕj   j   j 

31

(18)

The system of ODEs corresponding to the ‘iso-reactivity’ model is given in Table 2. The total reactivity of the protein is therefore the sum of the reactivities of all the residues, k0 = Nk and the effect of j , the doP, is predicted by the function ( 1 − j N ) ϕ j . Like this, the effect of the chemistry is decoupled from the diffusion/shielding effects arising from the PEG conjugation. The ‘iso-reactivity’ model can be useful in the case of more selective PEGylation strategies (i.e. often involving a single PEGylation site), as for cysteine PEGylation, to investigate the chemistry independently of the effect of the PEG molecular weight. Contrary to the ‘iso-reactivity’ model, the ‘infinite-selectivity’ model mixes these effects in apparent kinetic rate constants which are then losing their physical meaning. In the frame of the present model, even assuming ‘iso-reactivity’ of the lysine residues, the kinetic parameters retrieve their physical meaning and this is achieved with only three parameters, namely k , α and κ ' . Moreover, when looking at the data shown in Figure 8, this assumption seems reasonable. Indeed, for lysozyme a group of three residues with very close reactivities can be recognized, while the last three could be lumped together, so that in the end four sites of almost equal reactivity could be considered. Therefore, in the case of lysozyme PEGylation at the lysine residues, we obtained an experimental validation of the ‘iso-reactivity’ assumption. Table 2: Simplified kinetic model dP0 dt dPj dt

= −Nk ϕ0P0PEG

Disappearance of the native protein by PEGylation.

= (N − ( j − 1))k ϕ j −1Pj −1PEG − (N − j )k ϕ j Pj PEG

Appearance of the degree j and consumption to form the j +1.

dPEGd dPEG =− − dt dt

dPEGd dt

= kd PEG

N

∑j j =1

dPj dt

Disappearance of the activated PEG. Appearance of the hydrolysed PEG.

In Figure 11 is shown the comparison between simulations obtained with the ‘iso-reactivity’ model and experimental data. The simplified model was used to simulate the PEGylation of 32

lysozyme with three fold molar excess of PEG of various molecular weights. The results are presented in Figure 11 a-c. It can be seen that the model is in very good agreement with the experimental data, thus confirming the validity of the ‘iso-reactivity’ assumption. The best fitted values for this system are k =10.7 mL.mol-1.min-1, α =3.1 10-5 mol.g-1 and κ ' = 3.2 103 Å2. Moreover, literature data from Molek were also simulated (Molek, 2008). They correspond to the PEGylation of α-lactalbumin with 5-fold molar excess PEG5. The parameters for α-lactalbumin are k =1.26 mL.mol-1.min-1, α =9.0 10-5 mol.g-1 and κ ' = 0.0. The comparison between simulated and experimental values is shown in Figure 11d. It is worth noting that, in the present case where α-lactalbumin is PEGylated in only one condition, the ‘iso-reactivity’ model only requires two parameters, k and α , while the usual ‘infinite-selectivity’ model from Laura et al. would require five apparent kinetic parameters ( k1 =126, k2 =44, k3 =23, k4 =14 and k5 =9 mL.mol-1.min-1, from (Molek, 2008)). Therefore, the simplified ‘iso-reactivity’ model is proved to rely on reasonable assumptions leading to fairly good simulation results with a minimum number of fitting parameters. Moreover, the thus determined kinetic parameter, k is decoupled from the steric contributions ( α and κ ' ), so that the effect of the different contributions on the overall rate of PEGylation can be investigated independently.

33

Figure 11: Comparison between simulations and experimental data. The simulation were performed with the ‘iso-reactivity’ model. (a), (b) and (c) correspond to the PEGylation of HEWL (6 g/L solution) with 3-fold excess PEG5, PEG10 and PEG20, respectively. The best fitted values are k =10.7 mL.mol-1.min-1, α =3.1 10-5 mol.g-1 and κ ' = 3.2 103 Å2. (d) corresponds to the PEGylation of α-lactalbumin with 5-fold molar excess PEG5 (data from Molek (Molek, 2008)), with the parameters k =1.26 mL.mol-1.min-1, α =9.0 10-5 mol.g-1 and κ ' = 0.0. ((
) P0, () P1, () P2, () P3, () P4, () P5 and () PEG).

5.

Conclusion

In this work a kinetic model of protein PEGylation accounting for both PEGamers distribution and distribution of the positional isomers was developed. This provides a complete and fully detailed description of the product distribution resulting from the PEGylation process of a given protein. In addition, a specific model has been derived to compute the rate constant of the PEGylation reaction on the i th active site of a protein already conjugated with j PEG chains, k ji . The model 34

assumes that each i th active site on the protein exhibits its own reactivity which depends on the number j of PEG units already present in the protein but not on their position. The effect of the latter on the value of k ji is accounted for explicitly using an extended Smoluchowski model through two adjustable parameters: the shielding effect parameter α and the kinetic parameter κ ' . The conjugated PEG has been observed to wrap around the protein and form a ‘shell’ that

prevents some sites, at least temporarily, from being PEGylated. To describe this shielding effect, the kinetic rate constants were weighted by the fraction of accessible surface which was shown to scale logarithmically with the molecular weight of the PEG and the proportionality

(

constant α . On the other hand, the parameter κ ' 1 RPEG + 1 RP

j

) (R

c PEG

+ RPc

0

)

quantifies the

relative importance of the diffusion on the reaction rate. A large value of κ ' would imply that the chemical reaction of PEGylation is slow compared to the diffusion while a small value implies that the process is mainly controlled by the diffusion of the reacting species. The present model is general and can be used for any protein with N reactive sites. It requires in total N ‘intrinsic’ kinetic rate constants, two adjustable parameters, α and κ ' , and a hydrolysis rate constant for the functionalized PEG, kd . This PEG deactivation by hydrolysis of the ester bound is explicitly accounted for and described as a pseudo-first order kinetic reaction. A simplification of the kinetic model is proposed based on the assumption of equal reactivity of the PEGylation sites. This model requires only 3 kinetic parameters since the hydrolysis rate constant is estimated independently. From this point of view, the ‘iso-reactivity’ model allows to predict the effect of the doP, j with two parameters, α and κ ' , no matter the concentrations of reactants or the PEG molecular weight. The present models, simplified or fully detailed, can easily be used for rapid simulation or for optimizing the operating conditions while having a very good description of all the species involved when scaling up the process.

35

Acknowledgments This work was financially supported by Sanofi, Vitry sur Seine, France. In addition, we wish to thank our colleague Alberto Cingolani for his meticulous experimental work.

6.

References

Abuchowski, A., McCoy, J.R., Palczuk, N.C., van Es, T., Davis, F.F., 1977a. Effect of covalent attachment of polyethylene glycol on immunogenicity and circulating life of bovine liver catalase. J Biol Chem 252, 3582-3586. Abuchowski, A., van Es, T., Palczuk, N.C., Davis, F.F., 1977b. Alteration of immunological properties of bovine serum albumin by covalent attachment of polyethylene glycol. J Biol Chem 252, 3578-3581. Birch, J.R., Onakunle, Y., 2005. Biopharmaceutical Proteins. Humana Press. Calef, D.F., Deutch, J.M., 1983. Diffusion-Controlled Reactions. Annu Rev Phys Chem 34, 493524. Collins, F.C., Kimball, G.E., 1949. Diffusion-Controlled Reaction Rates. J Coll Sci Imp U Tok 4, 425-437. Day, J.N.E., Ingold, C.K., 1941. Mechanism and kinetics of carboxylic ester hydrolysis and carboxyl esterification. T Faraday Soc 37, 0686-0699. Einstein, A., 1905. Über die von der molekularkinetischen Theorie der Wärme geforderte Bewegung von in ruhenden Flüssigkeiten suspendierten Teilchen. Ann Phys 322, 549-560. Fee, C.J., Van Alstine, J.M., 2004. Prediction of the viscosity radius and the size exclusion chromatography behavior of PEGylated proteins. Bioconjugate Chem 15, 1304-1313. Fuchs, N., 1934. Über die Stabilität und Aufladung der Aerosole. Z Physik 89, 736-743. Hammes, G., 1978. Principles of chemical kinetics. . Academic Press. Harris, J.M., Kozlowski, A., 1997. Poly(ethylene glycol) and related polymers monosubstitued with propionic or butanoic acids and functional derivatives thereof for biotechnical applications, US005672662A ed. Shearwater Polymers, Inc., US. ICH, 2005a. Harmonized Tripartite Guideline: Pharmaceutical Development Q8. ICH, 2005b. Harmonized Tripartite Guideline: Quality Risk Management Q9. ICH, 2007. Draft Consensus Guideline Guideline: Pharmaceutical Development Annex to Q8. ICH, 2012. Harmonized Tripartite Guideline: Development and Manufacture of Drug Substances (Chemical Entities and Biotechnological/Biological Entities) Q11. Laura, B., Sa, H., Jianming, M., John, B., Rory, F., 2007. PEGylation of Biological Macromolecules, Process Chemistry in the Pharmaceutical Industry, Volume 2. CRC Press, pp. 383-402. Lee, B., Richards, F.M., 1971. The interpretation of protein structures: estimation of static accessibility. J Mol Biol 55, 379-&. Lee, H., Park, T.G., 2003. A novel method for identifying PEGylation sites of protein using biotinylated PEG derivatives. J Pharm Sci 92, 97-103. Maiser, B., Kroner, F., Dismer, F., Brenner-Weiss, G., Hubbuch, J., 2012. Isoform separation and binding site determination of mono-PEGylated lysozyme with pH gradient chromatography. J Chromatogr A 1268, 102-108. Meng, W., Guo, X., Qin, M., Pan, H., Cao, Y., Wang, W., 2012. Mechanistic Insights into the Stabilization of srcSH3 by PEGylation. Langmuir 28, 16133-16140. Molek, K.S., 2008. Ultrafiltration of PEGylated proteins. Pennsylvania State University. 36

Moosmann, A., Blath, J., Lindner, R., Muller, E., Bottingert, H., 2011. Aldehyde PEGylation Kinetics: A Standard Protein versus a Pharmaceutically Relevant Single Chain Variable Fragment. Bioconjugate Chem 22, 1545-1558. Moosmann, A., Christel, J., Boettinger, H., Mueller, E., 2010. Analytical and preparative separation of PEGylated lysozyme for the characterization of chromatography media. J Chromatogr A 1217, 209-215. Noyes, R.M., 1961. Effects of diffusion rates on chemical kinetics Prog React Kinet Mech 1, 129160. Noyes, R.M., 1969. Diffusion controlled reactions. Appl Spectrosc 23, 664-&. Pfister, D., Morbidelli, M., 2014. Process for protein PEGylation. J Control Release 180, 134149. Saltzman, W.M., 2001. Drug Delivery. Engineering Principles for Drug Therapy, Oxford University Press ed. Sandkühler, P., 2004. Kinetics of aggregation and gel formation in colloidal dispersions. ETH, Zürich. Schulten, K., 2000. Rate of Diffusion Controlled Reactions, Non-Equilibrium Statistical Mechanics. Smales, C.M., Moore, C.H., Blackwell, L.F., 1999. Characterization of lysozyme-estrone glucuronide conjugates. The effect of the coupling reagent on the substitution level and sites of acylation. Bioconjugate Chem 10, 693-700. Suckau, D., Mak, M., Przybylski, M., 1992. Protein surface topology-probing by selective chemical modification and mass-spectrometric peptide-mapping. Proc Natl Acad Sci USA 89, 5630-5634. Veronese, F.M., Mero, A., Pasut, G., 2009. Protein PEGylation, basic science and biological applications, in: Veronese, F. (Ed.), PEGylated Protein Drugs: Basic Science and Clinical Applications. Birkhäuser Basel, pp. 11-31. von Smoluchowski, M., 1906. Zur kinetischen Theorie der Brownschen Molekularbewegung und der Suspensionen. Ann Phys 326, 756-780. von Smoluchowski, M., 1917. Versuch einer mathematischen Theorie der Koagulationskinetik kolloidaler Lösungen. Z Phys Chemie 92, 129-168. Wang, X., Hu, J., Pan, D., Teng, H., Xiu, Z., 2014. PEGylation kinetics of recombinant hirudin and its application for the production of PEGylated HV2 species. Biochem Eng J 85, 38-48. Yang, C., Lu, D., Liu, Z., 2011. How PEGylation Enhances the Stability and Potency of Insulin: A Molecular Dynamics Simulation. Biochemistry 50, 2585-2593. Zalipsky, S., Harris, J.M., 1997. Introduction to chemistry and biological applications of poly(ethylene glycol). Amer Chemical Soc, Washington.

37

• • • • •

The amine PEGylation of lysozyme is investigated with various polymer sizes. A detailed kinetic model for protein PEGylation is described. Multiple reactive sites at the protein surface are accounted for. Diffusion-limited reactions and shielding effect are implemented. A limiting case assuming iso-reactivity of the sites is derived.

38