SSPBE: um programa para solução numérica da equação de Poisson-Boltzmann em simetria esférica com modelo de adsorção

SSPBE: a program for numerical solution of the spherically symmetric Poisson-Boltzmann equation along with a site binding model

Poisson-Boltzmann equation; site binding model; effective pair potential

Nota Técnica

SSPBE: UM PROGRAMA PARA SOLUÇÃO NUMÉRICA DA EQUAÇÃO DE POISSON-BOLTZMANN EM SIMETRIA ESFÉRICA COM MODELO DE ADSORÇÃO

Luís Gustavo Dias*, Hernan Chaimovich e Mário José Politi

Departamento de Bioquímica, Instituto de Química, Universidade de São Paulo, CP 26077, 05513-970 São Paulo - SP

*e-mail: [e-mail]

Recebido em 21/11/01; aceito em 3/5/02

SSPBE: A PROGRAM FOR NUMERICAL SOLUTION OF THE SPHERICALLY SYMMETRIC POISSON-BOLTZMANN EQUATION ALONG WITH A SITE BINDING MODEL. A Fortran77 program, SSPBE, designed to solve the spherically symmetric Poisson-Boltzmann equation using cell model for ionic macromolecular aggregates or macroions is presented. The program includes an adsorption model for ions at the aggregate surface. The working algorithm solves the Poisson-Boltzmann equation in the integral representation using the Picard iteration method.

Input parameters are introduced via an ASCII file, sspbe.txt. Output files yield the radial distances versus mean field potentials and average molar ion concentrations, the molar concentration of ions at the cell boundary, the self-consistent degree of ion adsorption from the surface and other related data.

Ion binding to ionic, zwitterionic and reverse micelles are presented as representative examples of the applications of the SSPBE program.

Keywords: Poisson-Boltzmann equation; site binding model; effective pair potential.

INTRODUÇÃO

O potencial de força média, medida de como as partículas interagem em um meio qualquer, é fundamental na compreensão de fenômenos de equilíbrio em sistemas iônicos e não-iônicos. Com o potencial de força média, as funções de distribuição e propriedades termoquímicas da solução podem ser calculadas¹. No caso de sistemas envolvendo macroíons (como micelas e vesículas carregadas, proteínas, polieletrólitos) e íons pequenos, como por exemplo, Na⁺, K⁺, Cl^-, Al⁺³, SO₄^-2, a obtenção do potencial de força média entre as partículas vai se tornando mais difícil à medida que a assimetria de carga e tamanho vai aumentando².

Nesta nota, é apresentada uma forma de se obter o potencial médio e a função de distribuição radial entre macroíons e microíons, através do programa "Spherically Symmetric Poisson-Boltzmann Equation (SSPBE)".

O programa SSPBE resolve numericamente a equação de Poisson-Boltzmann em simetria esférica no modelo de cela^3-5. Nesta abordagem, o volume da solução é dividido pelo número de macroíons, fornecendo o volume esférico médio ocupado por um macroíon na solução. Os contraíons e coíons são tratados como partículas puntiformes. O solvente é tratado como um dielétrico sem estrutura que blinda as interações eletrostáticas.

O programa permite trabalhar com três superfícies esféricas concêntricas carregadas de raios R₁, R₂ e R₃ (R₁£ R₂£ R₃) e os microíons podem ser distribuídos, segundo a estatística de Boltzmann, no volume existente entre as superfícies de raio R₁ e R₃ (Figura 1). A superfície esférica intermediária é permeável aos microíons e nenhum microíon pode existir fora da região definida pelo intervalo R₁£ r £ R₃. O programa também permite incluir um modelo de adsorção para microíons presentes nas superfícies esféricas.

Definido dessa forma, SSPBE é isomórfico a sistemas como uma solução que contém: a) micelas iônicas; b) micelas zwitteriônicas; c) micelas reversas com ou sem soluto no centro do compartimento aquoso; d) micelas mistas (iônica/zwitteriônica, ionica/não-iônica); e) vesículas impermeáveis a íons (negligenciando a distribuição iônica no espaço interior da bicamada quando tratamos o meio externo e vice-versa).

Na Figura 2, o esquema de uma micela reversa catiônica mais soluto catiônico hidrofílico é apresentado para exemplificar o isomorfismo. As três superfícies esféricas da Figura 1 estão representadas aqui, porém as superfícies 2 e 3 são coincidentes (R₂ = R₃), definindo o tamanho do compartimento aquoso da micela reversa. O soluto catiônico, tomado como esférico, aparece no centro da micela com raio R₁.

Todos estes sistemas podem ou não conter eletrólitos além dos contraíons provenientes do anfifílico.

O POTENCIAL DE CAMPO MÉDIO SEGUNDO A EQUAÇÃO DE POISSON-BOLTZMANN NA REPRESENTAÇÃO INTEGRAL

A lei de Coulomb⁶ aplicada ao sistema de três superfícies esféricas concêntricas resulta na equação 1 (em CGS):

onde: q(u) = densidade de carga total para microíons na dupla camada elétrica ("u" é a variável de integração e representa distância radial); q₁, q₂ e q₃ = carga das superfícies esféricas 1, 2 e 3; e = constante dielétrica efetiva (a mesma em todos os meios); R₁, R₂ e R₃ = raios das superfícies esféricas (R₁£ R₂£ R₃).

q(u) segue a distribuição de Boltzmann:

onde: r_i e = densidade numérica na superfície da cela e densidade numérica média da espécie i, respectivamente; Vc = volume da cela = 4p (R3 ¾ R1)³/3.

As densidades numéricas médias dos microíons são escolhidas de forma que a cela seja eletricamente neutra:

O MODELO DE ADSORÇÃO PARA ESTIMATIVA DO GRAU DE DISSOCIAÇÃO

Os microíons em solução podem ser adsorvidos nas superfícies esféricas. Ao ligarem nas superfícies, os microíons podem ser subdividos em fração ligada e fração livre. A adsorção de microíons pode ser estimada usando um modelo de ação de massas que relaciona o número de sítios na superfície esférica com a concentração de uma determinada espécie iônica⁷:

onde K_N,k = constante de equilíbrio de adsorção da espécie N na superfície k; n_k= número de sítios não ocupados na superfície k; n_Nk = número de sítios ocupados com o microíon da espécie N na superfície k; r_N(R_k)= densidade numérica da espécie N na superfície da casca k.

O modelo pode ser facilmente estendido para a adsorção de vários microíons com mais de uma superfície esférica, simplesmente usando o balanço para o número de sítios:

Também é possível definir a fração de carga da superfície k, a_k, resultante do balanço entre a quantidade de carga proveniente do número de microíons adsorvidos e a carga nominal da superfície:

onde z_k = valência do sítio localizado na superfície esférica k e z_i = valência da espécie i.

O balanço vale para cada superfície k, e n_total,k pode diferir de superfície para superfície. Além disso, é assumida a restrição de apenas um microíon por sítio.

A constante de equilíbrio pode ser estimada impondo-se a condição de concordância entre o valor para o grau de dissociação experimental (fração de carga da superfície) com o previsto pelo modelo na mesma situação.

Outros modelos de adsorção foram formulados similarmente ao tipo usado aqui^8,9.

SOLUÇÃO NUMÉRICA DA EQUAÇÃO DE POISSON-BOLTZMANN

A solução numérica passa pela discretização das integrais na equação 1 resultando em¹⁰:

onde:

_i = (r_i); q _j = q(r_j); r_i = r₁ + (N_in+1)h; N_in1 = número de intervalos usados na discretização das integrais entre as superfícies 1 e 2; N_in2 = número de intervalos usados na discretização das integrais entre as superfícies 2 e 3; h₁= espaçamento (uniforme) para integração entre as superfícies 1 e 2 = (R₂ - R₁)/N_in1; h₂ = espaçamento (uniforme) para integração entre as superfícies 2 e 3 = (R₃ - R₂)/N_in2; w_j = peso para a quadratura escolhida.

A discretização gera um sistema não-linear de equações para o potencial de campo médio. Uma alternativa para a solução da equação 7 é o uso do método iterativo simples (chamado também de método de Picard ou método de substituição sucessiva). Exemplificando para a equação 7:

O processo é iniciado escolhendo uma solução tentativa para o potencial médio, definido por

⁽⁰⁾, e a partir daí, calculando q_i⁽⁰⁾ da equação 2. Aplicando esta estimativa da densidade de carga no lado direito da equação 8, obtemos ⁽¹⁾. Para controle da convergência do processo iterativo, o potencial médio iterado (⁽ⁿ⁺¹⁾) é combinado com o potencial médio da iteração anterior (⁽ⁿ⁾) usando-se um parâmetro de mistura (g):

Isto resulta no potencial médio '⁽ⁿ⁺¹⁾ que será usado no processo iterativo. Em geral, g é próximo de zero e deve ser escolhido por tentativa e erro. Na nossa experiência, se g for muito pequeno, a convergência do processo iterativo é muito lenta e, se g for maior do que um certo valor crítico, o processo iterativo diverge. O processo é terminado quando o potencial médio usado numa iteração é próximo, segundo alguma definição de tolerância, ao potencial médio calculado.

A regra de quadratura implementada foi a trapezoidal. Nesta quadratura, os pesos são: w_i = 1, para i ¹ 1 e i ¹ N_in+1; w₁ = w_Nin+1 = ½.

O cálculo autoconsistente das frações de carga das superfícies esféricas com o modelo de adsorção também usa um método iterativo simples. A seqüência de procedimentos pode ser resumida em: a) definir inicialmente as frações de carga das superfícies; b) calcular iterativamente o potencial médio até convergência; c) usar o potencial médio nas equações 4, 5 e 6 para redefinir as frações de carga e, por conseqüência, as cargas nas superfícies esféricas; d) reiniciar o processo até a convergência das frações de carga.

Para evitar um processo divergente, é usado um parâmetro de mistura como aquele na equação 9. Algumas idéias recentes estendem este método simples de controle de convergência ao uso de um parâmetro de mistura que se ajusta de iteração em iteração¹¹.

COMPARAÇÃO COM CÁLCULOS DA LITERATURA E EXEMPLOS COM MICELAS

Em todos os sistemas discutidos abaixo, o tempo médio para cada cálculo em função da concentração do macroíon ou microíon foi de 1 à 2 min em um PC PentiumIII 500 MHz com sistema operacional Linux/GNU (Debian 2.2). Os cálculos foram efetuados em dupla precisão.

Micela iônica e vesícula carregada em solução de eletrólito 1:1

Mille e Vanderkooi³ como também Bell e Dunning⁴ resolveram numericamente a equação de Poisson-Boltzmann para esferas impermeáveis em simetria esférica, variando a concentração de sal e a concentração do macroíon. Outros grupos^12-15 fizeram estudos em simetria cilíndrica com o objetivo de comparar cálculos de atividade iônica e coeficiente osmótico com dados de macromoléculas helicoidais.

Os cálculos de Mille e Vanderkooi foram comparados com aqueles obtidos por SSPBE para dois sistemas: 1) micelas formadas de dodecilsulfato de sódio 0,05 M (CMC = 8 mM)³: as micelas iônicas são tratadas como esferas impermeáveis tendo um raio de aproximadamente 1,8 nm, número de agregação médio de 57, raio de cela igual a 8,2 nm; 2) vesículas formadas de ácido fosfatídico 1,9 mM³: as vesículas, também tratadas como esferas impermeáveis de raio igual a 30 nm, raio de cela igual a 172 nm e negligenciando a distribuição iônica no seu espaço interno, já que é considerado eletricamente neutro.

Para o modelo de solução micelar, as curvas de potencial reduzido, e/k_BT, em função da concentração de sal adicionado são mostradas na Figura 3. O potencial reduzido na superfície do macroíon e a razão entre as concentrações molares médias de contraíon/coíon são mostradas na Tabela 1, junto com os dados retirados de Mille e Vanderkooi.

Na Tabela 1 também estão os resultados para o modelo de cela da solução de vesículas. As curvas de potencial reduzido aparecem na Figura 4.

Os resultados apresentados na Tabela 1 indicam maior diferença entre os nossos cálculos e os de Mille e Vanderkooi para o caso da micela iônica (mesmo assim, não ultrapassam 0,6%).

Micela reversa variando o tamanho do compartimento aquoso

Karpe e Ruckenstein⁷ estudaram os efeitos de eletrólitos e tamanho do compartimento aquoso no potencial médio, concentração superficial de íons e grau de dissociação de micelas reversas de bis-(2-etilhexil)-sulfosuccinato de sódio, Aerossol OT, AOT, usando uma representação da equação de Poisson-Boltzmann na forma diferencial.

Os cálculos com SSPBE, abaixo apresentados, mostram que o potencial reduzido responde às mudanças nos parâmetros que definem as propriedades do compartimento aquoso (constante dielétrica, temperatura, tamanho, etc) de forma semelhante aos cálculos efetuados por Karpe e Ruckenstein⁷.

Nos cálculos variamos os tamanhos das micelas reversas simulando assim a dependência do tamanho da micela com a razão molar água-detergente (W/S). O número médio de agregação também é dependente de W/S. Os raios e números de agregação foram retirados de Karpe e Ruckenstein. O grau de dissociação foi mantido constante e igual a 0,28, provavelmente um limite superior ao valor real^7,16.

Os cálculos de potencial reduzido em função de W/S são apresentados na Figura 5. Também foram efetuados cálculos supondo a existência de uma molécula hidrofílica no centro da micela reversa (valência = -4, raio = 0,65 nm). Os valores simulam o ânion do 8-hidroxi-1, 3,6-trisulfonato-pireno (piranina). Este composto tem sido usado como sonda de atividade da água e de dinâmica de prótons em meios homogêneos e microheterogêneos^17,18.

Micela zwitteriônica

Embora as micelas zwitteriônicas sejam eletricamente neutras, possuem a capacidade de ligar e trocar íons (geralmente, ânions ligam mais fortemente do que cátions e o processo depende da orientação do grupo de átomos que definem o momento de dipolo permanente do "zwitterion")^19-21.

Em 1992, Baptista et al.²⁰ propuseram um modelo simples para tratar a interação de microíons com micelas zwitteriônicas. A micela foi pensada como um capacitor esférico com as superfícies representando a região dos dipolos da cabeça do detergente. A superfície mais externa é permeável aos microíons e água. Nesse caso, o modelo trata a ligação dos microíons como uma absorção pela micela.

SSPBE foi usado para estimar o efeito de sal no potencial reduzido de uma micela zwitteriônica de dodecildimetilamônio-propanosulfonato, DDPS. As micelas formadas por este anfifílico tem um número médio de agregação de 55 e concentração micelar crítica de 3,6 mM, com pequena alteração na presença de sal^18,22.

A partir do número de agregação, o raio da superfície esférica interna do capacitor é estimada em 1,8 nm e o da superfície externa em 2.4 nm. O efeito de sal no potencial reduzido em função da distância radial é apresentado na Figura 6 (linhas contínuas). O potencial reduzido da micela zwitteriônica em um meio sem sal aparece como uma linha pontilhada.

Na Figura 6 também estão curvas de potencial reduzido calculadas usando o modelo de adsorção como descrito anteriormente (linhas tracejadas). Nestes cálculos, o ânion do eletrólito interage de forma específica com a superfície interna positiva da micela (constante de dissociação, K, igual a 5 M). O potencial reduzido fica mais negativo em função do maior número de ânions presentes na região das cabeças do agregado, já que uma fração dos ânions encontra-se adsorvida na superfície positiva da micela.

CONCLUSÕES

O programa foi inicialmente desenvolvido para ser usado na estimativa de concentrações locais de contraíons e coíons, em interfaces micelares com uso em modelagem de reações termicamente ativadas²³ e obtenção do potencial médio de interação entre próton-fotoácido em micela reversa para a simulação do processo prototrópico da sonda (transientes de absorbância e fluorescência)¹⁸.

Embora limitado à simetria esférica, foi mostrado que SSPBE pode ser útil em uma série de problemas que envolvam colóides carregados e eletrólitos.

Os autores acreditam que SSPBE pode seu utilizado no ensino de graduação, em cursos que discutam teoria de eletrólitos.

O programa é gratuito e pode ser requisitado aos autores pelo correio eletrônico.

AGRADECIMENTOS

Os autores agradecem à FAPESP, CNPq e CAPES pelo apoio financeiro. Agradecem também à G. A. Marson, doutorando do IQUSP, pela ajuda com as Figuras 1 e 2 deste artigo. L.G. Dias é particularmente grato à FAPESP por uma bolsa de pós-doutoramento (Processo: 99/07688-9).

REFERÊNCIAS

1. Hansen, J.P.; McDonald, I.R.; Theory of Simple Liquids, Academic Press: London, 1986; McQuarrie, D.A.; Statistical Mechanics, Harper & Row: New York, 1976.

2. Lobaskin, V.; Lyubartsev, A.; Linse, P.; Phys. Rev. E 2001, 6302, 401; Hansen, J.P.; Lowen, H.; Annu. Rev. Phys. Chem. 2000, 51, 209; Vlachy, V.; Annu. Rev. Phys. Chem. 1999, 50, 145.

3. Mille, M.; Vanderkooi, G.; J. Colloid Interface Sci. 1977, 59, 211.

4. Bell, G. M.; Dunning, A. J.; Trans. Faraday Soc. 1970, 66, 500.

5. Marcus, R. A.; J. Chem. Phys. 1955, 23, 1057.

6. Reitz, J. R.; Milford, F. J.; Christy, R. W.; Foundations of Eletromagnetic Theory, Addison-Wesley: New York, 1980.

7. Karpe, P.; Ruckenstein, E.; J. Colloid Interface Sci. 1990, 137, 408.

8. Beunen, J. A.; Ruckenstein, E.; J. Colloid Interface Sci. 1983, 96, 469.

9. Feitosa, E.; Neto, A. A.; Chaimovich, H.; Langmuir 1993, 9, 702.

10. Kopchenova, N. V.; Maron, I. A.; Computational Mathematics, Mir: Moscou, 1987.

11. Busigin, A.; Phillips, C.R.; Mol. Phys. 1992, 76, 89.

12. Bell, G. M.; Dunning, A.; Trans. Faraday Soc. 1970, 66, 500.

13. Scheraga, H. A.; Katchalsky, A.; Alterman, Z.; J. Am. Chem. Soc. 1969, 91, 7242.

14. Rinaudo, M.; Loiseleur, B.; J. Chim. Phys. - Chim. Biol. 1973, 70, 1305.

15. Nitta, K.; Sugai, S.; J. Phys. Chem. 1974, 78, 1189.

16. Wong, M.; Thomas, J.K.; Nowak, T.; J. Am. Chem. Soc. 1977, 99, 4730.

17. Gutman, M.; Nachliel, E.; Annu. Rev. Phys. Chem. 1997, 48, 329.

18. Dias, L.G.; Tese de Doutorado, Universidade de São Paulo, Brasil, 1999.

19. Chevalier, Y.; Kamenka, N.; Chorro, M.; Zana, R.; Langmuir 1996, 12, 3225.

20. Baptista, M. S.; Cuccovia, I. M.; Chaimovich, H.; Politi, M. J.; Reed, W. F.; J. Phys. Chem. 1992, 96, 6442.

21. Baptista, M. S.; Politi, M. J.; J. Phys. Chem. 1991, 95, 5936; Brochsztain, S.; Berci-Filho, P.; Toscano, V. G.; Chaimovich, H.; Politi, M. J.; J. Phys. Chem. 1990, 94, 6781; Bunton, C. A.; Mhala, M. M.; Moffat, J.R.; J. Phys. Chem. 1989, 93, 854.

22. Herman, K. W.; J. Colloid Interface Sci. 1966, 22, 352.

23. Bunton, C. A.; Nome, F.; Quina, F. H.; Romsted, L. S.; Acc. Chem. Res. 1991, 24, 357.

Datas de Publicação

• Publicação nesta coleção
  20 Nov 2002
• Data do Fascículo
  Nov 2002

Histórico

• Aceito
  03 Maio 2002
• Recebido
  21 Nov 2001