Aplicação da equação de Poisson-Boltzmann ao cálculo de propriedades dependentes do pH em proteínas

Aplications of the Poisson-Boltzmann equation to the calculation of pH-dependent properties in proteins

non-linear Poisson-Boltzmann equation; pKa calculations; electrostatics in proteins

DIVULGAÇÃO

Aplicação da equação de Poisson-Boltzmann ao cálculo de propriedades dependentes do pH em proteínas

Aplications of the Poisson-Boltzmann equation to the calculation of pH-dependent properties in proteins

Thereza A. SoaresI, ^*; Ricardo Ferreira^II

^ILaboratory of Physical Chemistry, Swiss Federal Institute of Technology, ETH Hönggerberg, HCI G223, CH-8093, Zürich, Switzerland

^IIDepartamentos de Química Fundamental, Universidade Federal de Pernambuco, 50670-920 Recife - PE

ABSTRACT

The ability of biomolecules to catalyze chemical reactions is due chiefly to their sensitivity to variations of the pH in the surrounding environment. The reason for this is that they are made up of chemical groups whose ionization states are modulated by pH changes that are of the order of 0.4 units. The determination of the protonation states of such chemical groups as a function of conformation of the biomolecule and the pH of the environment can be useful in the elucidation of important biological processes from enzymatic catalysis to protein folding and molecular recognition. In the past 15 years, the theory of Poisson-Boltzmann has been successfully used to estimate the pK_a of ionizable sites in proteins yielding results, which may differ by 0.1 unit from the experimental values. In this study, we review the theory of Poisson-Boltzmann under the perspective of its application to the calculation of pK_a in proteins.

Keywords: non-linear Poisson-Boltzmann equation; pK_a calculations; electrostatics in proteins.

INTRODUÇÃO

O pH é uma propriedade físico-química essencial para a viabilidade de sistemas biológicos, e inúmeras doenças estão associadas com variações do pH sanguíneo em apenas 0,5 unidade do valor fisiológico 7,4. Os efeitos nocivos resultantes destas variações não são causados pela ação química de íons H₃O⁺ ou OH^- no sangue. Se assim fosse, um pH sanguíneo de 6,8 não seria fatal visto que uma variação tão pequena de pH não resulta em troca significante de prótons entre biomoléculas e ambiente. Uma explicação mais correta para este fato seria a influência do pH sobre a estabilidade e reatividade de biomoléculas, dado que alterações do pH correspondentes a uma energia térmica RT são fatais para seres humanos. Considerando-se que DDG = DpH(RT/log e) = DpH(RT/0,43), uma variação do pH em 0,4 unidades é crítica para sistemas biológicos¹. E de fato, este valor define perfeitamente o intervalo no qual alterações do pH sanguíneo são toleradas por seres humanos.

As moléculas biológicas são particularmente sensíveis a variações de pH porque são constituídas por diversos grupos químicos que ionizam em pH aproximadamente igual a 7,0. Por exemplo, uma dada proteína pode muito bem ser estável em pH 5,0 ou 9,0 desde que histidinas, cisteínas e grupos aminos N-terminais estejam ausentes. O motivo pelo qual a Seleção Natural escolheu biomoléculas tão sensíveis a variações de pH está, em parte, relacionado com a importância de reações de transferência de prótons em catálise enzimática. A condição para que um grupo químico participe eficientemente em uma reação de transferência de prótons é que o mesmo seja apto a doar ou receber prótons. Adicionalmente, este grupo químico deve ser capaz de retornar prontamente ao estado de ionização anterior à transferência de prótons, regenerando assim o estado enzimático ativo. Para que um grupo químico preencha estes requisitos, o mesmo deve possuir um pK_a aproximadamente igual ao pH do ambiente. Conseqüentemente, catálise enzimática é necessariamente sensível às variações do pH.

A dependência entre a estrutura e a função de biomoléculas com relação ao pH é um desafio importante à compreensão e manipulação de propriedades biomoleculares. Em proteínas, a variação de pH é de origem predominantemente eletrostática². Assim, a determinação de cargas fracionárias de grupos químicos em função de propriedades como conformação e pH do ambiente seria de enorme valor para a elucidação de inúmeros processos biológicos. Como exemplo podemos citar, além dos mecanismos de catálise enzimática, o enovelamento de proteínas e o reconhecimento molecular. Diversos resíduos constituintes da estrutura nativa de proteínas apresentam um comportamento de ionização diferente daquele quando isolados em solução, devido à influência do ambiente eletrostático da proteína sobre o grupo ionizável.

A grande maioria dos aminoácidos que constituem uma proteína interage através de forças de longo alcance, essencialmente de natureza eletrostática. Os campos elétricos gerados por proteínas estendem-se a uma distância de 10-15 Å, dependendo de condições de temperatura, solvente e carga elétrica da proteína (a maiores distâncias, as flutuações térmicas aleatórias do solvente predominam sobre qualquer efeito que o campo elétrico da proteína possa exercer). Desta forma, para entender as forças de atração-repulsão entre biomoléculas é necessário compreender as leis eletrostáticas que regem as interações intermoleculares. Os fundamentos da física eletrostática são bem estabelecidos e podem ser expressos concisamente por equações relativamente simples. Contudo, a aparente simplicidade destas equações esconde algumas dificuldades conceituais e numéricas da aplicação das mesmas ao estudo de sistemas complexos. Este é um problema relevante devido à vasta quantidade de estruturas tridimensionais de biomoléculas disponíveis recentemente.

Este artigo apresenta uma introdução ao uso da equação de Poisson-Boltzmann na determinação de estados de protonação em proteínas. Antes de prosseguir, é importante lembrar que como qualquer modelo teórico, o modelo de Poisson-Boltzmann é uma simplificação da realidade e envolve o uso de várias aproximações. Assim, a aplicação da equação de Poisson-Boltzmann ao calculo pK_a em proteínas pressupõe as seguintes premissas como verdadeiras: i) as contribuições eletrostáticas e não-eletrostáticas em um dado processo químico podem ser calculadas separadamente, ii) princípios macroscópicos podem ser aplicados ao nível molecular. Aqui, assume-se que a diferença entre o comportamento de protonação de um dado grupo químico em solução, para o qual a constante de ionização é conhecida, e do mesmo grupo químico no ambiente protéico, é de natureza puramente eletrostática. O cálculo das respectivas energias livres eletrostáticas é obtido através da solução da equação de Poisson-Boltzmann pelo uso do método das diferenças finitas³. O enfoque do artigo é limitado aos modelos que tratam o soluto e o solvente como meios dielétricos contínuos. Desta forma, métodos baseados em simulações microscópicas, os quais são importantes para o estudo de moléculas orgânicas pequenas, não são abordados nesta revisão, mas uma literatura extensiva está disponível na área⁴.

TEORIA

A equação de Poisson-Boltzmann

O cálculo de energias eletrostáticas entre grupos químicos em uma proteína pode ser obtido através da equação de Poisson-Boltzmann linearizada. Para isto é necessária a determinação dos potenciais eletrostáticos associados com uma determinada distribuição de cargas. Uma descrição clássica destes potenciais é dada pela equação de Poisson, que mostra como varia a densidade de carga r da atmosfera iônica com a distância do íon na posição r. Se a distribuição de cargas tivesse simetria esférica, a equação de Poisson tomaria a forma:

No caso mais geral, a equação de Poisson pode ser escrita:

A equação de Poisson permite relacionar a variação espacial do potencial f na posição r com a distribuição da densidade de carga r em um meio de constante dielétrica e. O operador Ñ (nabla) no sistema de coordenadas cartesianas tem a forma (¶/¶x, ¶/¶y, ¶/¶z). No contexto de um sistema biológico composto por cargas imersas em meio aquoso contendo eletrólitos, a carga parcial de cada átomo na molécula pode ser descrita como cargas fixas r^f, enquanto que as cargas dos eletrólitos dissolvidos no solvente podem ser descritas como cargas móveis, cuja densidade é determinada pela distribuição de Boltzmann:

onde n_i é a densidade de íons do tipo i por unidade de volume em uma dada região do espaço; n⁰_i a densidade de íons do tipo i por unidade de volume em todo o espaço amostrado; q_i a carga correspondente a cada íon; f o potencial eletrostático em uma dada região do espaço; k_B a constante de Boltzmann e T a temperatura absoluta. Assim, a Equação (1) pode ser modificada na chamada equação de Poisson-Boltzmannn:

onde l é igual a 1 em regiões acessíveis aos íons móveis e 0 nas demais regiões. O efeito da força iônica do meio é tratado através da teoria de Debye-Hückel, a qual permite a aplicação da equação de Poisson-Boltzmannn a sistemas que podem ser adequadamente representados como um conjunto de cargas em um meio de constante dielétrica uniforme e de baixa força iônica. Em tais condições é possível tratar linearmente a distribuição de cargas móveis representada pelo termo exponencial na Equação (3):

onde k é a constante de Debye-Hückel ou blindagem iônica do meio; e a carga eletrônica; N_a o número de Avogadro e I a força iônica. As variáveis F, e, k e r são funções do vetor posição r. A forma da interface dielétrica e a exclusão de íons contidos no interior da proteína são introduzidos no cálculo através da dependência das constantes e e k do vetor r. Em condições onde k = 0, íons estão ausentes e a Equação (4) pode ser reduzida à equação de Poisson. Devido à diferença de polarizabilidade da água e do interior da proteína, duas constantes dielétricas - e e e_p - são requeridas para representar realistamente o sistema.

Propriedades dielétricas

O tratamento de biomoléculas em meio contínuo é baseado nos modelos desenvolvidos por Born⁵, Onsager⁶ e Tanford e Kirkwood⁷, nos quais o soluto é descrito como uma cavidade de baixa constante dielétrica e imerso em um meio de constante dielétrica e=78. Nesta aproximação, a proteína e o composto modelo são representados como um objeto rígido de constante dielétrica e_p=4-6. Ambas moléculas são compostas por sítios i, cuja protonação-deprotonação é representada pela adição-subtração de uma carga unitária pontual ("single-site model"), e por sítios j, que representam cargas permanentes na molécula. A interface entre soluto-solvente é definida pelo contorno resultante do movimento de um átomo de raio definido ao longo da superfície de van der Waals dos átomos constituintes da proteína ou do composto modelo.

Um aspecto fundamental para o uso apropriado de modelos de meio contínuo é a representação adequada de diferentes ambientes dielétricos através da definição de constantes dielétricas. Para isto, a resposta do ambiente ao campo elétrico gerado pela proteína ou pelo composto modelo deve ser considerada. Esta resposta pode ser dividida em três fenômenos físicos que resultam em um escalonamento de cargas atômicas: i) redistribuição de cargas, ii) polarizabilidade eletrônica e iii) reorientação de dipolos permanentes (moléculas polares).

Quanto a redistribuição de cargas na aproximação de campo-médio, os íons não são tratados explicitamente. Ao invés disto, assume-se que o potencial químico de cada íon seja uniforme ao longo da solução. Em um meio de constante dielétrica uniforme e baixa força iônica, a redistribuição de cargas pode ser descrita apropriadamente pelo termo de Debye-Hückel na equação de Poisson-Boltzmannn (Equação 4). Um tratamento mais sofisticado da mobilidade iônica ao redor da proteína ou do composto modelo é possível através do uso de simulações do tipo Monte Carlo^8-10, onde efeitos de correlação e do tamanho finito dos íons são considerados. Contudo, resultados similares podem ser obtidos com os dois métodos, Monte Carlo e Poisson-Boltzmann, para a distribuição de íons univalentes e de tamanho reduzido¹¹. Nos casos de íons divalentes, a equação de Poisson-Boltzman apresenta um desempenho quantitativo mais fraco, embora qualitativamente tais resultados permanecem válidos¹¹.

O segundo efeito, a polarizabilidade eletrônica é o resultado da reorientação da nuvem eletrônica ao redor do núcleo atômico, quando em presença de um campo elétrico. Em líquidos orgânicos, este processo freqüentemente resulta em constantes dielétricas de 1,5-2. Devido à energia eletrostática e constante dielétrica serem grandezas inversamente proporcionais, e fazendo e=2, um soluto imerso em um líquido orgânico apolar apresenta uma energia eletrostática que corresponde à metade da energia do mesmo soluto no vácuo. Por este ser um problema de múltiplos-corpos, a modelagem do mesmo é complexa e computacionalmente dispendiosa. Como conseqüência, a maioria dos campos de força empregados em simulações do tipo mecânica e dinâmica molecular não incluem explicitamente os efeitos resultantes de polarização eletrônica¹². Várias alternativas para o tratamento de polarizabilidade eletrônica em meios contínuos estão disponíveis, e estes modelos foram revisados extensivamente por Sharp e Honig¹³. O modelo dielétrico uniforme é comumente acoplado ao uso da equação de Poisson-Boltzmann, e os efeitos resultantes de polarizabilidade eletrônica podem ser incorporados ao valor da constante dielétrica, definindo um dado ambiente eletrostático. Este modelo assume que cargas e dipolos constituindo o sistema interagem entre si como se estivessem imersos em ambiente de constante dielétrica 2. Uma questão relevante é a validade do uso de um parâmetro único e invariável para descrever uma região constituída por inúmeros e diversos tipos de átomos. Os argumentos em favor do modelo dielétrico uniforme - comumente usado pelo modelo de Poisson-Boltzmann - são baseados em resultados experimentais. Primeiramente, a constante dielétrica de alta freqüência característica de líquidos orgânicos varia muito pouco com relação à identidade da molécula de solvente. Em segundo lugar, modelos mais sofisticados, nos quais constantes dielétricas de diferentes grupos químicos seriam considerados, requerem a disponibilidade de valores experimentais de polarizabilidades para um número largo de átomos em diferentes moléculas. Por outro lado, o modelo dielétrico uniforme requer apenas o conhecimento da quantidade experimental e (disponível para inúmeros solventes orgânicos) para representar a polarizabilidade eletrônica em sistemas tão complexos como proteínas.

As propriedades dielétricas de líquidos polares resultam da re-orientação de dipolos permanentes em um campo elétrico. Para um sistema composto por uma proteína imersa em meio aquoso é necessário considerar os dipolos permanentes correspondentes ao soluto e ao solvente. Assim, a simulação destes sistemas exige a representação dos dipolos permanentes da macromolécula e do solvente através do uso de duas constantes dielétricas. O estudo experimental de propriedades dielétricas de biomoléculas é extremamente difícil e, conseqüentemente, a validação de diferentes métodos teóricos torna-se inacessível¹³. Diferentes valores de constante dielétrica são utilizados em estudos teóricos de proteínas. Estes valores são obtidos através da aplicação direta de equações, tais como a equação de Coulomb^14,15, o modelo de Tanford-Kirkwood^7,16, a equação de Poisson-Boltzmann^3,17, ou através de teoria dielétrica de Kirkwood-Frölich^18-21. Estes estudos sugerem que a constante dielétrica de proteínas e_p possuem valores médio em torno de 4-6.

Idealmente, constantes dielétricas são utilizadas como forma de reproduzir aspectos do sistema que não são tratados explicitamente. No caso de modelos contínuos, e_p sumariza as contribuições devidas a efeitos de polarizabilidade eletrônica, reorientação de dipolos e aspectos estruturais que não são inclusos explicitamente na representação estrutural da proteína ou interface proteína-solvente. Estes aspectos estruturais incluem rearranjos relacionados com a redistribuição de cargas no interior e na superfície da proteína, bem como flutuações do sistema em um dado estado de carga. A comparação entre diferentes valores de e_p indica que a escolha de e_p=20 resulta em valores de pK_a mais próximos aos valores determinados experimentalmente^22-24. Numa situação ideal, rearranjos estruturais estariam acoplados a mudanças no campo elétrico de modo que interações entre grupos carregados são otimizadas; uma constante dielétrica alta compensa a ausência de flexibilidade estrutural inerente ao modelo. Assim, na ausência de modelos mais detalhados, a constante dielétrica pode ser considerada como um parâmetro ajustável e definido através da comparação entre resultados teóricos e experimentais^25,26.

Desenvolvimento do modelo eletrostático para o cálculo de pK_a em proteínas: uma perspectiva histórica

A primeira aplicação do modelo de meio contínuo ao cálculo de pK_as em biomoléculas foi realizada em 1924 por Lindestrøm-Lang²⁷. O chamado "smeared-charge model" trata a proteína como um íon de Born gigantesco, com a carga líquida de todos os grupos ionizáveis espalhada ao longo da superfície da esfera. Em 1957, Tanford e Kirkwood⁷ introduziram modificações ao modelo através do tratamento dos grupos ionizáveis como cargas discretas em posições arbitrariamente escolhidas ao longo da superfície da proteína representada por uma esfera de baixa constante dielétrica imersa em solvente, representado por um meio contínuo de alta constante dielétrica. O efeito da força iônica sobre o sistema foi incorporado ao modelo, permitindo a redistribuição dos eletrólitos dissolvidos no meio em resposta ao potencial eletrostático da proteína, de acordo com a equação linear de Poisson-Boltzmannn. No entanto, a principal limitação destes modelos foi a inexistência de estruturas atômicas de proteínas, juntamente com a capacidade computacional limitada da época. À medida que informações sobre a estrutura cristalográfica de proteínas e valores experimentais de pK_as para aminoácidos em solução tornaram-se disponíveis, as equações do modelo de Tanford-Kirkwood foram utilizadas para estimar interações entre grupos ionizáveis, incluindo os valores cristalográficos recém obtidos para a distância entre pares de átomos constituindo grupos ionizáveis²⁸.

O modelo de Tanford-Kirkwood permaneceu limitado pela representação simplista de proteínas como esferas e pela aproximação de que o pK_a de um dado grupo ionizável na proteína é perturbado unicamente por interações com outros grupos ionizáveis. Conseqüentemente, a influência de grupos neutros, dipolos e o efeito de desolvatação são negligenciados.

Em 1982, Warwicker e Watson³ propuseram um modelo matemático para solucionar a equação de Poisson-Boltzmann aplicada a um sistema proteína-solvente, empregando o método das diferenças finitas. Este modelo possibilitou tratar a forma da proteína mais precisamente, a qual passou a ser representada pela superfície de van der Waals, e introduziu um termo representativo das interações entre o grupo ionizável e grupos neutros e/ou dipolos associados com a presença de motivos de estrutura secundária. Contudo, este modelo continuou desprezando a energia de desolvatação, isto é, a contribuição resultante do trabalho eletrostático de transferir um grupo polar de um ambiente solvatado - representado por uma constante dielétrica alta - para o interior da proteína - representado por uma constante dielétrica baixa.

Em 1988, este problema foi resolvido através do cálculo do potencial eletrostático f_i para cada átomo i carregado ou parcialmente carregado em um dado grupo ionizável no ambiente do solvente e no interior da proteína¹⁸. A energia eletrostática total ou energia de desolvatação é dada pelo trabalho resultante da transferência do grupo ionizável entre os dois ambientes dielétricos.

Métodos numéricos

Diversos métodos numéricos foram desenvolvidos para solucionar a equação de Poisson-Boltzmannn aplicada a sistemas biológicos^2,29,30. O método das diferenças finitas é o modelo matemático implementado em programas como DELPHI, UHBD e MEAD entre outros^{13,19,29,31-35}. Nesta aproximação, a macromolécula e os seus arredores correspondendo ao solvente são representados por uma caixa reticulada, onde cada ponto localizado no vértice do retículo corresponde a uma porção da proteína ou do solvente. Valores de densidade de carga, constante dielétrica, força iônica e potencial eletrostático são assinalados à cada ponto do retículo. A equação de Poisson-Boltzmannn linearizada é, então, integrada para todo o volume da caixa reticulada³. Na forma mais simples do modelo, a reação de ionização é representada pela adição de uma carga ±1 a um átomo que recebe ou perde um próton durante a reação ("single-site model"). Este átomo corresponde ao sítio de ionização. Por exemplo, os átomos nitrogênio N^z e carbono C^g são os sítios de ionização em uma lisina e um ácido aspártico, respectivamente. A equação de Poisson-Boltzmannn linearizada é então resolvida para a forma protonada e deprotonada de cada sítio de ionização no sistema. Primeiramente, a equação é resolvida para o grupo imerso no ambiente dielétrico da proteína, a qual por sua vez está imersa no ambiente dielétrico do solvente. O limite entre a proteína e o solvente é tipicamente a superfície molecular ou de van der Waals da proteína. Este primeiro cálculo permite avaliar o campo de reação devido à influência do solvente sobre a carga do sítio ionizável, bem como a interação do mesmo com outras cargas na proteína. Subseqüentemente, o grupo ionizável é transferido do ambiente dielétrico da proteína para aquele do solvente, possibilitando a determinação do campo de reação produzido pelo solvente quando o grupo está completamente solvatado. O modelo assume que o pK_a do grupo quando imerso no solvente é igual ao valor do pK_a determinado experimentalmente para os aminoácidos em solução^36-37.

Termodinâmica do equilíbrio de protonação

O equilíbrio de protonação do grupo químico AH isolado em solução AH_s e presente em uma proteína AH_p pode ser representado por um ciclo termodinâmico conforme descrito na Figura 1.

A variação da energia livre padrão de deprotonação do grupo ionizável AH em solvente DGº_s é dada pela equação:

onde DGº_s é igual à soma entre a energia livre de deprotonação do grupo AH em fase gasosa DGº_g, as energias livres de solvatação do grupo na forma protonada DG_s(AH) e deprotonada DG_s(A^-) e a energia livre de solvatação do próton DG_s(H⁺):

Assim, o pK_a do grupo AH em um solvente pode ser determinado a partir da variação da energia livre padrão de deprotonação DGº_s quando os valores para os demais termos são conhecidos. Os valores das energias livres padrão em fase gasosa podem ser obtidos através de métodos ab-initio^4,38 ou medidas experimentais³⁵ (e.g. espectroscopia de ressonância magnética ou cinética enzimática), enquanto que as energias livres de solvatação podem ser calculadas por métodos Monte Carlo^39,40 ou perturbação termodinâmica⁵⁶. Contudo, a energia livre de solvatação do próton DG_g,s(H⁺) é uma quantidade difícil de ser determinada, tanto experimental quanto computacionalmente, de modo que a relação resultante da diferença entre os pK_as de dois grupos AH e BH é preferivelmente expressa como:

Analogamente, o pK_a do grupo AH na proteína, pK_a^(p), é descrito pela equação:

Os valores de pK_a^(s) são obtidos a partir de medidas experimentais para cada um dos 20 aminoácidos naturais em solução³⁵. A Equação (9) pode ser rearranjada na seguinte equação:

A Equação (10) é utilizada preferencialmente na implementação computacional porque, de modo similar aos cálculos de energia livre, a determinação de diferenças de energia entre estruturas semelhantes é mais simples e precisa.

A relação entre o pK_a, o pH e a constante de equilíbrio K_a para a reação de deprotonação de um sítio ionizável é dada pela equação de Henderson-Hasselbach:

O rearranjo das Equações (6) e (11) permite calcular a variação da energia livre em função do pH e determinar a curva de titulação para um sítio ionizável:

Teoria do equilíbrio de protonação

Estritamente, reações de ionização de grupos químicos não são processos puramente eletrostáticos, visto que dependem implicitamente do potencial químico do próton e da energia de ligação do próton aos grupos ionizáveis. Contudo, efeitos de pH são largamente mediados por alterações de estado de ionização, isto é, alteração de cargas elétricas em grupos químicos ácidos e básicos. Assim, interações eletrostáticas podem ser consideradas efetivamente como forças primárias controlando processos químicos dependentes do pH^2,7. Conseqüentemente, a diferença entre a energia livre de protonação de um grupo ionizável em uma proteína e a energia livre de protonação do mesmo grupo em solução pode ser descrito apropriadamente pela diferença entre o trabalho eletrostático de alterar a carga do grupo no estado desprotonado para o estado protonado na proteína e o trabalho eletrostático para a mesma alteração no composto modelo^22,32,41,42.

Assumindo-se que pares de grupos ionizáveis no sistema interagem de forma aditiva, é possível representar as energias livres de interação entre estes grupos por uma matrix simétrica ||G_ij||. Considerando o estado de referência para cada grupo ionizável como aquele onde o grupo está na forma neutra, podemos definir os termos da matrix de energia livre ||G_ij||. Cada termo diagonal (i=j) representa a energia livre resultante da diferença entre o trabalho de deprotonar o grupo i quando todos os demais grupos ionizáveis estão na forma neutra e o trabalho de deprotonar o mesmo grupo em solução. Os termos fora da diagonal (i¹j) representam as energias livres de interação quando os grupos i e j estão simultaneamente desprotonados. Assim, a diferença entre a energia livre de deprotonação de um sistema com M sítios ionizáveis (s₁, ,s_M ) pode ser representado em termos da matrix de energia livre ||G_ij||:

onde x_i é igual a 1 quando o grupo está na forma desprotonada e igual a 0 quando o grupo está na forma neutra, g_i é igual a +1 para bases e -1 para ácidos.

Neste modelo, a diferença entre a energia eletrostática de ionização de um grupo na proteína e do mesmo grupo em solução é determinada por três tipos de contribuições (Figura 2): i) a energia de desolvatação de Born, DDG_Born, é a energia resultante da interação entre o grupo ionizável e a polarização induzida nos ambientes dielétricos da proteína e do solvente, respectivamente, pela carga deste grupo; ii) a energia de interação de "background", DDG_Back, descreve a interação entre o grupo ionizável e as cargas parciais presentes na proteína, quando todos os outros grupos ionizáveis estão na forma neutra. Este termo inclui também a interação resultante da presença de dipolos associados com motivos estruturais como a-hélices na vizinhança do grupo ionizável. Estas duas contribuições correspondem aos termos G_ii na matrix de energias ||G_ij|| como mostrado na Equação (13); iii) a energia de interação do grupo ionizável com os demais grupos ionizáveis na proteína. Este termo é calculado exclusivamente para o grupo ionizável no ambiente dielétrico da proteína, sendo representado pelos termos G_ij na matrix de energias ||G_ij||.

Contudo, o cálculo das energias de interação entre um dado grupo ionizável e os demais grupos ionizáveis na proteína é complexo, porque as cargas destes grupos são, por si mesmas, dependentes do pH. A introdução de um termo hipotético pK_i permite contornar esta dificuldade⁷. Este termo corresponde ao pK_a do grupo ionizável na proteína, quando os demais grupos ionizáveis estão na forma neutra. A diferença entre o pK_a de um grupo ionizável e o pK_i da proteína pode ser determinada pelo cálculo de DDG_Born e DDG_Back:⁷

Os elementos da matrix de energia livre de interação ||G_ij|| podem ser calculados pelas Equações (15) e (16):

onde f_i é o potencial eletrostático ao redor da carga q_i. Este potencial é a soma entre o campo de reação devido à carga q_i e todas as interações coulômbicas resultantes da presença da carga q_i entre as demais cargas presentes nos dois ambientes dielétricos em consideração (solvente e proteína). Na Equação (14), o termo q_j representa a carga do grupo j, f é o potencial resultante da interação entre as cargas dos grupos j e i. Os estados protonados são marcados por um asterisco (*), em oposição aos estados neutros do mesmo grupo ionizável. Os primeiros dois termos correspondem aos potenciais resultantes da interação entre a forma carregada do grupo i e as cargas das formas carregadas e neutras do grupo j, respectivamente. Os dois últimos termos correspondem aos potenciais resultantes da interação entre a forma neutra do grupo i e as cargas das formas carregadas e neutras do grupo j, respectivamente. Os termos presentes nas Equações (13) e (14) representam o somatório dos produtos entre cargas e potenciais na vizinhança destas cargas e, conseqüentemente, representam energias livres eletrostáticas.

Tratamento de sistemas com múltiplos sítios de protonação

A modelagem do equilíbrio de protonação em proteínas é particularmente difícil devido ao grande número de sítios de ionização e de interações eletrostáticas de longo-alcance entre estes sítios. Em proteínas, cada um dos M grupos ionizáveis pode ocorrer em pelo menos dois estados de protonação diferentes (ou mais para grupos com diferentes tautômeros, como a histidina). Em outras palavras, a determinação de propriedades dependentes do pH em uma proteína requer o cálculo das energias eletrostáticas para todos os possíveis estados de protonação do sistema e subseqüente avaliação dos valores médios destas propriedades baseada na distribuição de Boltzmann. No entanto, o número de possíveis estados de protonação é muito alto para a maioria dos sistemas biológicos. É necessário o uso de aproximações que possibilitem a determinação de valores médios para uma dada propriedade, sem a enumeração de todos os possíveis estados de protonação do sistema em consideração.

Diferentes métodos foram propostos para tratar o problema de protonação de múltiplos sítios^{1,32,41,43,44}. Entre estes, o método de Monte Carlo baseado no algoritmo de Metropolis é o mais comumente utilizado. Uma proteína apresenta um total de 2^M estados de ionização possíveis, onde M é o número de sítios ionizáveis presentes no sistema. Esta simplificação implica que cada sítio de ionização apresenta apenas dois estados possíveis: protonado ou desprotonado^37,38,41.

O estado de ionização inicial corresponde àquele onde todos os sítios ionizáveis estão nas respectivas formas protonadas. A energia livre G1 para o estado de protonação inicial pode ser calculada através da equação:

A Equação (23) é uma expressão modificada da Equação (19), na qual o estado de referência para cada grupo ionizável corresponde à forma protonada dos mesmos. Na Equação (19), o estado de referência para cada grupo ionizável corresponde à forma neutra destes grupos. A variável x_i' é igual a 1 quando o grupo está na forma protonada e 0 quando está na forma deprotonada e q_i⁰ é a carga formal do sítio desprotonado. Subseqüentemente, um número no intervalo entre 0 e 1 é aleatoriamente assinalado à cada um dos sítios ionizáveis do sistema, e o estado de protonação é então alterado para todos os grupos ionizáveis cujo número randômico corresponde a um valor maior que um dado valor s. Assim, um número médio de (1-s)/100% dos grupos ionizáveis terão os respectivos estados de protonação alterados. Se s é igual a 0,99 aproximadamente 1% dos estados de protonação serão alterados em cada passo da simulação de Monte Carlo. O segundo estado de protonação é então aceito ou rejeitado dependendo da energia livre G2, novamente calculada pela Equação (23); quando G2 < G1, o novo estado de protonação é aceito, quando G2 > G1, o novo estado de protonação é aceito de acordo com a probabilidade baseada no valor de e^-(G2-G1)/RT. Um número aleatório no intervalo entre 0 e 1 é novamente escolhido, e se o mesmo é menor que e^-(G2-G1)/RT, o novo estado de protonação é aceito. A aplicação do método de Monte Carlo ao cálculo de múltiplos sítios de protonação é computacionalmente caro, mas é o método de escolha para tratamento de sistemas com numerosos átomos e cujos grupos ionizáveis sofrem deprotonação/protonação de modo acoplado¹.

Alguns fatores limitantes da precisão do modelo de Poisson-Boltzmann e possíveis alternativas para o problema

Apesar dos fundamentos teóricos do modelo de Poisson-Boltzmann terem sido definidos, vários fatores introduzem imprecisões significativas e limitam a capacidade de previsão do modelo³⁵. Como exemplo, pode-se citar a parametrização de cargas atômicas e raios de van der Waals (este último determinando a superfície molecular da proteína, i. e., a interface proteína-solvente), a representação do ambiente dielétrico da proteína como um meio homogêneo, o tratamento de moléculas de água posicionadas no interior da proteína e, principalmente, o tratamento da flexibilidade conformacional de grupos ionizáveis.

Inúmeras soluções foram propostas na tentativa de anular o efeito destes fatores sobre a capacidade de predição do método. Stikoff et al.⁴⁵ parametrizaram cargas e raios atômicos para reproduzir as energias de solvatação de pequenas moléculas orgânicas. Brucoleri et al.⁴⁶ introduziram uma aproximação na qual cargas atômicas são tratadas como esferas de carga uniforme, ao invés de cargas pontuais. A vantagem desta aproximação é menor dependência do método das diferenças finitas com relação ao posicionamento das cargas na caixa reticulada. Oberoi e Allewell⁴⁸ obtiveram melhor concordância entre resultados teóricos e calculados, ao assinalar a carga resultante da protonação/deprotonação ao átomo mais exposto ao solvente. Outra aproximação promissora é a determinação de raios atômicos com base na distribuição radial de moléculas de água ao redor dos respectivos átomos⁴⁹. Yang et al.⁴¹ mostraram que a inclusão de moléculas de água na vizinhança de grupos ionizáveis resulta em uma representação mais realística das propriedades de titulação em lisozima. Gibas e Subramaniam⁵⁰ determinaram sistematicamente o número de moléculas de água que deveriam ser inclusas explicitamente durante o cálculo para reproduzir apropriadamente os valores experimentais de pK_a em lisozima. Estes autores concluíram que aproximadamente uma molécula de água deve ser inclusa dentro do raio de uma ponte de hidrogênio para cada grupo ionizável na proteína.

Antosiewicz et al.²² propuseram o uso de uma constante dielétrica interna e_p=20, muito maior do que o valor e_p=2-4, comumente aceito para proteínas. Em princípio, uma constante dielétrica alta compensaria a omissão de variações conformacionais locais e da ligação específica de íons, conforme discutido na seção Propriedades dielétricas. Independentemente, Schaefer et al.²⁴ investigaram o uso de constantes dielétricas e_p variando de 1 a 30, e observaram que e_p = 20 resulta em valores de pK_a mais próximos ao valores experimentais. States e Karplus⁵¹ propuseram o uso de diferentes valores de e_p para representar resíduos no interior e na superfície da proteína. Este modelo foi posteriormente testado por Demchuck e Wade²³, resultando em melhor concordância entre os valores calculados e experimentais de pK_a. A vantagem deste modelo é a possibilidade de um tratamento fisicamente mais realístico de ambientes dielétricos tão distintos, quanto o interior e a superfície de uma proteína. Mais recentemente, foi proposto que na ausência de modelos mais detalhados para representar a permitividade interna da proteína, a constante dielétrica deveria ser considerada como um parâmetro ajustável, cujo valor deve ser definido através da comparação entre resultados teóricos e experimentais e de acordo com o sistema em consideração^25,26.

Indubitavelmente, a restrição mais grave da teoria de Poisson-Boltzmann para o cálculo de pK_a em biomoléculas é a ausência de flexibilidade conformacional do sistema. Pequenos rearranjos estruturais de grupos ionizáveis no sistema podem causar modificações significativas do ambiente eletrostático na vizinhança do sítio de ionização e, conseqüentemente, dos valores de pH estimados. Parte do problema é devido à própria inadequação de estruturas cristalográficas em representar um "ensemble" de conformações, o qual vem sendo parcialmente aliviado pelo número crescente de estruturas de proteínas obtidas por técnicas de ressonância magnética nuclear (RMN). Contudo, é importante ressaltar que estruturas determinadas por espectroscopia de RMN representam, de fato, uma conformação média e não uma única conformação, visto que a posição de cada átomo em tais estruturas não pode ser inferida diretamente como no caso de métodos cristalográficos. Além desta limitação, o refinamento das diferentes soluções espaciais que podem corresponder ao conjunto de intensidades espectroscópicas ("NOE restraints") é obtido por métodos de Mecânica Molecular, o qual faz uso dos chamados campos de força ("force-fields") parametrizados para biomoléculas em pH fisiológico, i.e., valores de pH entre 7,0-7,2.

Inúmeros métodos computacionais propuseram a inclusão de flexibilidade conformacional do sistema no cálculo de pK_a, pelo uso de um "ensemble" de estruturas obtidas por métodos como Monte Carlo ou Dinâmica Molecular^33,47,52-54. Em particular, cadeias laterais localizadas no interior da proteína representam um grande desafio ao modelo, porque freqüentemente apresentam valores de pK_a que diferem em várias unidades dos respectivos valores em solução. Em adição, conformações diferentes de uma mesma cadeia lateral resultam em valores calculados de pK_a que diferem grandemente. Uma solução para o problema é o cálculo do valor médio de pK_a para diferentes conformações do sistema, obtido por métodos de Dinâmica Molecular ou Monte Carlo^26,35,55. No entanto, as conformações obtidas por Dinâmica Molecular ou Monte Carlo são, por sua vez, dependentes dos estados de protonação escolhidos para a conformação inicial do sistema. Recentemente, métodos de dinâmica molecular em pH constante foram propostos, nos quais os estados de protonação podem variar ao longo da simulação^56,57. Estes métodos tratam o pH como um parâmetro termodinâmico, de modo similar à pressão e temperatura ("proton bath"), sem fazer uso da equação de Poisson-Boltzmann e estão fora do âmbito deste artigo. No entanto, estes são métodos muito caros computacionalmente e, até o momento, restritos a sistemas de tamanho reduzido, como é o caso do acido succínico.

CONCLUSÃO

A teoria de Poisson-Boltzmann permite o tratamento de propriedades dependentes do pH em proteínas, com um erro médio de 1 unidade de pK_a. A grande vantagem desta aproximação é o baixo custo computacional, comparativamente a outros métodos como integração termodinâmica, por exemplo. Em adição, a aplicação extensiva deste modelo a inúmeros sistemas permitiu uma caracterização aprofundada das limitações inerentes ao método. No entanto, é importante observar-se que a capacidade de previsão e a precisão do modelo de Poisson-Boltzmann dependem intimamente de fatores externos ao método, como a qualidade dos valores experimentais utilizados tanto durante o cálculo de pK_a quanto durante a validação dos valores calculados contra valores experimentais. Infelizmente, a grande maioria das proteínas possuem apenas um pequeno número de grupos ionizáveis cujos valores de pK_a são conhecidos experimentalmente e, freqüentemente, estes valores estão restritos às proteínas de tamanho reduzido, devido às limitações da técnica de RMN, comumente usada na determinação experimental de valores de pK_a em proteínas.

Do ponto de vista teórico, o modelo de Poisson-Boltzmann pode ser aperfeiçoado através da inclusão de alguns aspectos até então ignorados, como: i) modelos mais detalhados para representar a redistribuição de cargas, devido ao processo de protonação/deprotonação; ii) inclusão de flexibilidade ao menos das cadeias laterais pertencentes aos grupos ionizáveis; iii) inclusão de moléculas de água localizadas no interior da proteína e uso de uma constante dielétrica apropriada para representar tais moléculas e iv) representação apropriada da heterogeneidade dielétrica do ambiente protéico. Certamente, como em qualquer modelo teórico, o modelo de Poisson-Boltzmann é uma aproximação matemática que procura representar sistemas fisicamente complexos. As deficiências discutidas anteriormente podem ser melhoradas ao preço de maior custo computacional e desenvolvimento de novos algoritmos computacionais; uma das metas na área é o acoplamento entre o cálculo de estados de ionização e dinâmica molecular em proteínas. No entanto, como tem sido demonstrado pela extensa literatura na área, o modelo de Poisson-Boltzmann é uma importante ferramenta computacional no cálculo de pK_a em proteínas e tem sido amplamente utilizado para o tratamento realístico de propriedades dependentes do pH nestes sistemas.

REFERÊNCIAS

1. Gilson, M. K.; Proteins 1993, 15, 266.

2. Daves, M. E.; McCammon, J. A.; Chem. Rev. 1990, 95, 509.

3. Warwicker, J.; Watson, H. C.; J. Mol. Biol. 1982, 157, 671.

4. Evanseck, J. D.; Houk, K. N.; Briggs, J. M.; Jorgensen, W. L.; JACS 1994, 116, 10630; Jorgensen, W. L.; Briggs, J. M.; JACS 1989, 111, 4190; Jorgensen, W. L.; Briggs, J. M.; Gao, J.; JACS 1987, 109, 6857.

5. Born, M.; Z. Phys Chem. 1920, 1, 45.

6. Onsager, L.; JACS 1936, 58, 1486.

7. Tanford, C.; Kirkwood, J.; JACS 1957, 79, 5333; Tanford, C.; JACS 1957, 79, 5340.

8. Lebret, M.; Zimm, B. H.; Biopolymers 1984, 23, 271; Lebret, M.; Zimm, B. H.; Biopolymers 1984, 23, 287.

9. Mezei, M.; Beveridge, D. L.; Method. Enzymol. 1986, 127, 21.

10. Anderson, C. F.; Record, M. T.; Annu. Rev. Biophys. Biomol. Struct. 1990, 19, 423.

11. Murthy, C. S.; Bacquet, R. J.; Rossky, P. J.; J. Phys. Chem. 1985, 89, 701.

12. Gilson, M. K.; Curr. Opin. Struct. Biol. 1995, 5, 216.

13. Sharp, K. A.; Honig, B.; Annu. Rev. Biophys. Chem. 1990, 19, 301; Sharp, K. A.; Honig, B.; J. Phys. Chem. 1990, 54, 301.

14. Hill, T. L.; J. Phys. Chem. 1956, 60, 253.

15. Rees, D. C.; J. Mol. Biol. 1980, 141, 323.

16. Matthew, J. B.; Annu. Rev. Biophys. Biomol. Struct. 1985, 14, 387.

17. Honig, B.; Nicholls, A.; Science 1995, 268, 1144.

18. Gilson, M. K.; Honig, B.; Biopolymers 1986, 25, 2097.

19. Simonson, T.; Perahia, D.; Comput. Phys. Commun. 1995, 91, 291.

20. King, G.; Lee, F. S.; Warshel, A.; J. Chem. Phys. 1991, 95, 4366.

21. Smith, P. E.; Brunne, R. M.; Mark, A. E.; van Gunsteren, W.F.; J. Phys. Chem. 1993, 97, 2009.

22. Antosiewicz, J.; McCammon, J. A.; Gilson, M. K.; J. Mol. Biol. 1994, 238, 415.

23. Demchuk, E.; Wade, R. C.; J. Phys. Chem. 1996, 100, 17373.

24. Schaefer, M.; Sommer, M.; Karplus, M.; J. Phys. Chem. B 1997, 101, 1663.

25. Baptista, A. M.; Soares, C. M.; J. Phys. Chem. B 2001, 105, 293.

26. Soares, T. A.; Lins, R. D.; Straatsma, T. P.; Briggs, J. M.; Biopolymers 2002, 65, 313; Soares, T. A.; Goodsell, D. S.; Ferreira, R.; Olson, A. J.; Briggs, J. M.; J. Mol. Recognit. 2000, 13, 146.

27. Lindestrøm-Lang, K.; Comptes Rendus des Travaux du Laboratoire Carlsberg 1924, 15, 1.

28. Imoto, T.; Biophys. J. 1983, 44, 293.

29. Honig, B.; Sharp, K.; Yang, A.; J. Phys. Chem. 1993, 97, 1101.

30. Zhou; H. X.; J. Chem. Phys. 1994, 100, 3152.

31. Gilson, M. K.; Sharp, K.A.; Honig, B.; J. Comput. Chem. 1987, 9, 327.

32. Bashford, D.; Karplus, M.; Biochemistry 1990, 29, 10219; Bashford, D.; Karplus, M.; J. Phys. Chem. 1991, 95, 9556.

33. Bashford, D.; Gerwert, K.; J. Mol. Biol. 1992, 224, 473.

34. Gilson, M. K. Em Computer simulation of biomolécular systems; van Gusteren, W. F.; Weiner, P. K.; Wilkinson, A. J., eds.; Kluwer Academic Publishers: The Netherlands, 1997, vol. 3, p. 119-222.

35. Briggs, J. M.; Antosiewicz, J.; Em Reviews in Computational Chemistry; Lipkowitz, K. N.; Boyd, D. B., eds.; John Wiley & Sons: New York, 1999, vol. 13, p. 249-311.

36. Daves, M. E.; Madura, J. D.; Luty, B. A.; McCammon, J. A.; Comput. Phys. Commun. 1991, 62, 187.

37. Madura, J. D.; Briggs, J. M.; Wade, R. C.; Daves, M. E.; Luty, B. A.; Ilin, A.; McCammon, J. A.; Comput. Phys. Commun. 1995, 91, 57.

38. Custódio, R.; Politi, J. R. D.; Segala, M.; Haiduke, R. L. A.; Cyrillo, M.; Quim. Nova 2002, 25, 159; Freitas, L. C. G.; Quim. Nova 1999, 22, 293.

39. Barlette, V. E.; Freitas, L. C. G.; Quim. Nova 1999, 22, 254.

40. Freitas, L. C. G.; Silva, L. B.; Botelho, L. D.; Quim. Nova 1996, 19, 166.

41. Yang, A.; Gunner, M. R.; Sampogna, R.; Sharp, K.; Honig, B.; Proteins 1993, 15, 252.

42. Del Buono, G. S.; Figueirido, F. E.; Levy, R. M.; Proteins 1994, 20, 85.

43. Tanford, C.; Roxby, R.; Biochemistry 1972, 11, 2192.

44. Antosiewicz, J.; Porschke, D.; Biochemistry 1989, 28, 10077.

45. Sitkoff, D.; Lockhart, D. J.; Sharp, K. A.; Honig, B.; Biophys. J. 1990, 67, 2251.

46. Bruccoleri, R. E.; J. Comput. Chem. 1993, 14, 1417.

47. Beroza, P.; Fredkin, D. R.; Okamura, M. Y.; Feher, G.; PNAS 1991, 88, 5804.

48. Oberoi, H.; Allewell, N. M.; Biophys. J. 1993, 65, 48.

49. Nina, M.; Beglov D.; Roûx, B.; J. Phys. Chem. B 1997, 101, 5239.

50. Gibas, C. J.; Subramaniam, S.; Biophys. J. 1996, 71, 138.

51. States, D. J.; Karplus, M.; J. Mol. Biol. 1987, 197, 122.

52. Yang, A. S.; Honig, B.; J. Mol. Biol. 1993, 231, 459.

53. Bashford, D.; Case, D. A.; Dalvit, C.; Tennant, L.; Wright, P. E.; Biochemistry 1993, 32, 8045.

54. Antosiewicz, J.; Briggs, J. M.; Elcock, A. H.; Gilson, M. K.; McCammon, J. A.; J. Comput. Chem. 1996, 17, 1633.

55. Langsetmo, K.; Fuchs, J. A.; Woodward, C.; Sharp, K. A.; Biochemistry 1991, 30, 7609.

56. Børjesson, U.; Hünenberger, P. H.; J. Chem. Phys. 2001, 114, 9706.

57. Baptista, A. M.; Teixeira, V. H.; Soares, C. M.; J. Chem. Phys. 2002, 117, 4184.

Recebido em 6/5/03; aceito em 8/10/03; publicado na web em 27/05/04

Datas de Publicação

• Publicação nesta coleção
  31 Jul 2004
• Data do Fascículo
  Ago 2004

Histórico

• Recebido
  06 Maio 2003
• Aceito
  08 Out 2003