IMPORTANCIA DA MODELAGEM MOLECULAR PARA MOLECULAS BIOATIVAS

SUMÁRIO

1 INTRODUÇÃO................................................................................................11

2 MODELAGEM MOLECULAR........................................................................12

2.1 Processos modernos no desenvolvimento de fármacos.............................13

2.2 Métodos empregados na elaboração dos programas de modelagem molecular...........................................................................................................15

2.2.1 Mecânica molecular..................................................................................16

2.2.2 Dinâmica Molecular..................................................................................17

2.2.3 Métodos ab initio.......................................................................................18

2.2.4 Deformação no comprimento de ligação..................................................19

2.2.5 Barreira de energia de rotações intramoleculares (ângulos de torsão)................................................................................................................20

2.2.6 Interações de van der Waals....................................................................20

2.2.7 Métodos semi-empíricos...........................................................................20

2.2.8 Funcional de Densidade...........................................................................21

2.2.9 Métodos comparativos..............................................................................22

3 APLICAÇÃO DE MODELAGEM MOLECULAR EM QUÍMICA MEDICINAL.......................................................................................................22

3.1 Estudos de Inibidores de Cicloxigenase (Cox) e 5-Lipoxigenase (5-LO): Proposta de novos inibidores Seletivos de 5- LO..............................................23

3.2 Estudo Teórico da Tromboxana A2: Proposta de uma Conformação Bioativa..............................................................................................................29

4 MÉTODOS DE MODELAGEM PARA ESTUDO E PLANEJAMENTO DE COMPOSTOS BIOATIVOS...............................................................................36

5 CONCLUSÃO.................................................................................................38

REFERÊNCIAS BIBLIOGRÁFICAS.................................................................39

1 INTRODUÇÃO

Os altos custos de tempo e dinheiro são fatores preocupantes para o crescimento da Indústria Farmacêutica, que necessita de investimentos em Pesquisa e Desenvolvimento (P&D) para a descoberta de novos fármacos com maior eficácia e seletividade. Estima-se um gasto entre 800 milhões a 1,4 bilhões de dólares e cerca de 15-25 anos para o desenvolvimento de um novo fármaco (GELDENHUYS et. al, 2006). Assim, para otimizar tempo e custo, a utilização de inovações científicas e tecnológicas, como ferramenta de pesquisa, combinando conhecimentos multidisciplinares de informática, biotecnologia, química e biologia, refletem no planejamento de novos fármacos (DREWS, 2000; DREWS, 2003).

Planejamento de fármacos auxiliado por computador (Computer-Aided Drug Design, CADD) é uma dessas evoluções tecnológicas promissoras no desenvolvimento de novos candidatos a fármacos (TANG et. al., 2006), podendo reduzir em até 50% os gastos na pesquisa de uma novo agente terapêutico (GELDENHUYS et. al., 2006). Há muitas técnicas de CADD que auxiliam na descoberta de novas moléculas candidatas a fármacos e suas estratégias de utilização dependem das informações estruturais disponíveis do alvo (enzima/receptor) e do ligante (MEEK et. al, 2006).

A concepção de novos fármacos por CADD baseia-se em duas estratégias principais: planejamento “direto” e “indireto”. A primeira estratégia considera as características tridimensionais de um alvo conhecido, e.g. enzima ou receptor, analisando-se o complexo formado entre o ligante e a macromolécula-alvo, identificando-se possíveis sítios preferenciais de ligação e principais modos de interação. Estes dados permitem propor modificações específicas na estrutura do ligante, com o objetivo de aumentar a afinidade e especificidade ao alvo (VAN GUNSTEREN et. al., 1994; AJAY & MURCKO, 1995).

A segunda estratégia é usada quando a estrutura da macromolécula-alvo não é conhecida. Assim, informações sobre a atividade e características estruturais e estéreoeletrônicas dos compostos ativos e inativos podem ser utilizadas para determinar propriedades específicas de uma molécula que podem influenciar na interação com o alvo, tais como calor de formação, potencial eletrostático molecular (Molecular Eletrostatic Potential - MEP), densidade eletrônica, energia e coeficiente dos orbitais de fronteira HOMO (Highest Occupied Molecular Orbital) e do LUMO (Lowest Unoccupied Molecular Orbital), energia de ionização, ordem de ligação e momento de dipolo (BARREIRO et. al., 1997), grupos hidrofóbicos, grupos aceptores e doadores de ligação hidrogênio. A partir dessas informações, gera-se um modelo que pode ser utilizado para a seleção de compostos de bancos de dados ou orientar o processo de planejamento (COHEN, 1990).

A modelagem molecular reúne um conjunto de técnicas computacionais que possibilitam a construção, a visualização, a manipulação e a estocagem de modelos moleculares tridimensionais. A modelagem molecular permite a análise conformacional, o cálculo de propriedades estéreo-eletrônicas e a análise de variações estruturais que auxiliam na interpretação das correlações entre as estruturas químicas de uma série de compostos com a variação da atividade farmacológica, sendo de grande importância no planejamento de fármacos.

2 MODELAGEM MOLECULAR

A Modelagem Molecular envolve cálculos de parâmetros selecionados para determinada molécula, fornecendo informações sobre suas propriedades (Cohen, 1996).

Os métodos de estimativa ou cálculo de propriedades moleculares podem ser divididos em (Cohen, 1996):

• Interpolação: a correlação é encontrada entre a propriedade ou característica desejada e outras propriedades ou características moleculares;

• Extrapolação: a correlação não inclui a molécula de interesse, mas extrapolando a correlação é possível estimar a propriedade desejada.

• Cálculos computacionais: cálculos das propriedades moleculares realizadas por computador, com níveis diversos de variáveis, que podem substituir as medidas experimentais quando o procedimento necessário é muito elaborado ou de difícil realização (Cohen, 1996).

Algumas propriedades como: energia de formação, energia de ativação, momento de dipolo, características de orbitais moleculares, eletrônicas, atômicas, de ligações, polarizabilidade, caráter hidrofóbico, e outras como coeficiente de partição (logP) e refratividade molar (MR) podem ser originados por através de cálculos computacionais que incluem mecânica molecular, métodos semiempíricos, cálculos ab initio e outros métodos empíricos (Cohen, 1996; Sant’Anna, 2002).

Os três principais métodos teóricos para cálculos computacionais de propriedades podem ser classificados em métodos empíricos (mecânica molecular), métodos semi-empíricos e cálculos ab initio (Cohen, 1996).

Mecânica molecular e métodos semi-empíricos envolvem cálculos por meio de parametrizações e cálculos ab initio reproduzem um experimento sem utilizar parâmetros.

Os melhores resultados obtidos nos dois primeiros métodos são aplicados em interpolações com dados experimentais (Cohen, 1996; Leach, 2001).

Cálculos semi-empíricos são baseados em mecânica quântica, envolvendo funções de onda que incluem orbitais de Slater (STO) e de Gaussian (GTO). Métodos semi-empíricos tratam linearmente combinações dos orbitais interagindo computacionalmente com um campo de força estável (SCF) e permitindo a minimização de energia do sistema (Cohen, 1996).

Os métodos semi-empíricos diferem entre si pelas aproximações ao considerar a repulsão dos elétrons de diferentes orbitais. Estas aproximações são realizadas com valores de parametrização correspondentes aos obtidos por cálculos ab initio e/ou experimentais. Os métodos recentes usam parametrizações oriundas de cálculos ab initio baseados na teoria de Huckel estendida e no programa CNDO (Complete Neglet of Differential Overlap), que auxiliam nos cálculos de orbitais e cargas. A fim de obter dados mais seguros quanto à energia, métodos semiempíricos aperfeiçoados foram introduzidos como MINDO-3 (Modified Intermediate Neglect of Differential Overlap), MNDO (Modified Neglect of Diatomic Overlap) e AM1 (Austin Model 1) (Cohen, 1996; Sant’Anna, 2002)

Baseando-se nestas propriedades e em estudos prévios que correlacionam à estrutura química com a atividade biológica e o pesquisador em modelagem molecular analisa a série de moléculas para as quais efetuou os cálculos computacionais e pode predizer a atividade esperada. As análises experimentais serão necessárias para a validação da hipótese aventada.

2.1 Processos modernos no desenvolvimento de fármacos

Um dos mais importantes avanços no planejamento e descoberta de novos fármacos tem sido a utilização da modelagem molecular. Ela tem se firmado como uma ferramenta indispensável não somente no processo de descoberta de novos fármacos, mas também na otimização de um protótipo já existente ou obtido pelo próprio estudo de modelagem molecular.

O grande desenvolvimento da modelagem molecular deveu-se em grande parte ao avanço dos recursos computacionais em termos de hardware (velocidade de cálculo) e software (programas de modelagem molecular). No passado, a utilização desta técnica era restrita a um seleto grupo de pessoas que desenvolviam os seus próprios programas de modelagem molecular. Atualmente, não é mais necessário ao usuário (modelista) compor o seu próprio programa em virtude deles poderem ser obtidos através de grandes companhias e de laboratórios acadêmicos.

A modelagem molecular fornece informações importantes para o processo de descoberta de fármacos. Ela permite a obtenção de propriedades específicas de uma molécula que podem influenciar na interação com o receptor. Como exemplos, podemos citar o mapa de potencial eletrostático, o contorno da densidade eletrônica e a energia e os coeficientes dos orbitais de fronteira HOMO (Highest Occupied Molecular Orbital) e do LUMO (Lowest Unoccupied Molecular Orbital) etc. Outras informações importantes também podem ser obtidas a partir da comparação estrutural entre diferentes moléculas, o que pode permitir a geração de um índice de similaridade que pode ser correlacionado com a atividade farmacológica. A modelagem molecular também permite a visualização tridimensional (3D) do complexo fármaco-receptor e fornece informações sobre os requisitos estruturais essenciais que permitem uma interação adequada do fármaco no seu sítio receptor. Esta ferramenta também tem o potencial de planejar teoricamente novas moléculas que satisfaçam as propriedades eletrônicas e estruturais para um perfeito encaixe no sítio receptor.

A maioria dos programas de modelagem molecular é capaz de desenhar a estrutura molecular e realizar os cálculos de otimização geométrica e estudos de análise conformacional. Os arquivos de saída destes cálculos podem ser utilizados como arquivos de entrada para outros programas. Desta forma, a primeira etapa em estudos de modelagem molecular é desenhar a estrutura da molécula. Em seguida, a molécula é otimizada objetivando encontrar parâmetros geométricos tais como comprimentos e ângulos de ligação (Figura 1) que estejam próximos aos valores determinados experimentalmente. Desta forma, pode-se avaliar a qualidade do programa de modelagem molecular selecionado para efetuar os cálculos considerando que ele deve ser capaz de representar corretamente a estrutura molecular sem que os parâmetros estruturais da referida molécula tenham sido usados para elaborá-lo

Assim, um programa de modelagem molecular deve ser capaz de adotar o princípio da transferibilidade, ou seja, reconhecer e transferir os parâmetros embutidos no programa para uma nova molécula que apresente as mesmas características estruturais e eletrônicas das moléculas usadas para confeccionar o programa (mesmo tipo de átomos, funções químicas, hibridização molecular etc).

Figura 1: Representação do comprimento de ligação (d21 e d32) e ângulo de ligação (a321) definidos pelos átomos At1, At2 e At3.

A rotação entre ligações saturadas, como entre átomos de carbono, permite que uma única molécula adote diversas conformações. O estudo de análise conformacional permite determinar as conformações de mínimo de energia (confôrmeros). Estes confôrmeros indicam como os grupamentos funcionais estão orientados e, portanto, revelam aspectos relevantes de como a molécula pode interagir com um receptor específico, considerando que a conformação mais estável deve estar em maior número durante o processo de interação com o receptor. Entretanto, devemos ressaltar que não existe uma relação entre a conformação mais estável e a conformação bioativa, pois a primeira pode sofrer mudanças na sua conformação original no momento em que se aproxima do sítio receptor. Este estudo conformacional é realizado considerando as rotações livres entre dois átomos consecutivos que, por sua vez, determinam um ângulo diedro, assinalando como devem estar orientados os grupamentos adjacentes a estes dois átomos (Figura 2).

Figura 2: Representação do ângulo de diedro delimitados por quatro átomos consecutivos, assinalando diferentes valores de ângulo de diedro.

2.2 Métodos empregados na elaboração dos programas de modelagem molecular

Um programa de modelagem molecular permite a representação, visualização, manipulação e determinação de parâmetros geométricos (comprimento e ângulo de ligação) e eletrônicos (energia dos orbitais de fronteira, momento de dipolo, potencial de ionização etc) de uma molécula isolada, além de realizar estudos em macromoléculas (proteínas) e complexos droga–receptora.

A grande maioria dos programas de modelagem molecular é capaz de retratar entidades químicas com um alto grau de precisão. Esta afirmação é oriunda de estudos comparativos de parâmetros eletrônicos e geométricos obtidos experimentalmente. Os métodos mais empregados para a obtenção de propriedades moleculares são: a mecânica molecular e os métodos semi-empíricos.

2.2.1 Mecânica molecular

Mecânica molecular (MM) é um método que calcula a estrutura e a energia das moléculas com base nos movimentos dos núcleos. Os elétrons não são considerados explicitamente, mas, ao contrário, é assumido que eles encontrarão uma distribuição ótima, uma vez que as posições dos núcleos são conhecidas. Esta idéia é baseada na aproximação de Born-Oppenheimer. Esta aproximação estabelece que os núcleos são mais pesados e, portanto, movem-se mais lentamente do que os elétrons. Desta forma, os movimentos nucleares, as vibrações e as rotações podem ser estudadas separadamente,admitindo que os elétrons movem-se rapidamente e ajustam-se aos movimentos do núcleo. Assim, pode-se admitir que a mecânica molecular trata a molécula como uma coleção de esferas conectadas por molas, onde as esferas representam os núcleos e as molas representam as ligações (Figura 3).

O campo de força é usado para calcular a energia e a geometria de uma molécula. Ele é elaborado de forma que contenha uma coleção de diferentes tipos de átomos, parâmetros (para comprimento e ângulos de ligação e etc.) e equações para calcular a energia de uma molécula. Em um determinado campo de força, um dado elemento pode ter diferentes tipos de átomos. Por exemplo, o etilbenzeno contém átomos de carbono com hibridização sp3 e átomos de carbono aromático (sp2). Os átomos de carbono com hibridização sp3 apresentam uma geometria de ligação tetraédrica, enquanto átomos de carbono aromático têm uma geometria trigonal (planar).

Figura 3: Representação de uma molécula utilizando princípios de modelagem molecular onde as esferas são os átomos e amola representa a ligação entre eles.

O comprimento de ligação entre os átomos de carbono (C-C) no grupo etil difere do comprimento C-C no grupo fenila e o comprimento da ligação entre o átomo de carbono etila e o átomo de carbono da fenila ligados diretamente é diferente de todos os comprimentos de ligação no etilbenzeno. Os campos de força contêm parâmetros para estes diferentes tipos de ligação. A energia total de uma molécula pode ser dividida em várias partes denominadas forças potenciais ou equações de energia potencial. Estas forças potenciais são calculadas independentemente e somadas para obter a energia total da molécula. Exemplos de forças potenciais são as equações de energias associadas com a deformação do comprimento (s), ângulo de ligação (a), ângulo de torsão (t), interações de van der Waals (vdW) etc. Estas equações definem a superfície de energia potencial de uma molécula.

ETOTAL = Es + Ea + Et + Evdw

2.2.2 Dinâmica Molecular

Em sistemas reais, a energia cinética dos átomos permite que colisões intermoleculares e que mudanças conformacionais aconteçam a todo momento. Um modelo de trabalho que considera apenas a conformação de menor energia encontrada em um processo de minimização de energia de uma molécula bioativa tem limitações, visto que na temperatura da biofase a energia cinética disponível permitiria que muitas outras conformações fossem facilmente alcançadas.

A resolução das equações de movimento da dinâmica molecular representa a evolução no tempo dos movimentos moleculares, chamada trajetória, que pode ser usada para o estudo de propriedades que dependem do tempo, tais como a difusão, o dobramento de cadeias moleculares(como acontece no enovelamento de proteínas) e a distribuição de moléculas de solvente ao redor de um soluto, entre outros. Em um cálculo de dinâmica molecular, as moléculas, que têm uma certa energia cinética, podem superar barreiras de energia potencial, o que permite a exploração da superfície de energia potencial (SEP) na busca de outras conformações (ou distribuições, no caso de conjuntos de moléculas) estáveis.

Os cálculos de dinâmica molecular baseiam‐se na resolução da equação do movimento de Newton para cada átomo i do sistema molecular, onde Fi é a força que causa a aceleração al em um átomo de massa mi. A força pode ser obtida diretamente da derivada da energia potencial E com relação às coordenadas ri dos átomos

O tratamento clássico só permite a resolução destas equações para sistemas com até duas partículas independentes. A resolução para sistemas maiores é feita como uso de métodos numéricos. Se a posição no tempo t é r(t), a posição após um intervalo de tempo curto Δt pode ser obtida pela seguinte série de Taylor:

A solução numérica então depende do conhecimento da posição r(t), da velocidade dr/dt e da aceleração d²r/dt² para cada átomo, além de aproximações para as contribuições de termos de maior ordem. O intervalo de tempo Δt de cada etapa da trajetória é um parâmetro de grande importância em um cálculo de dinâmica molecular. Primeiramente, Δt deve ser pequeno o suficiente para que a aceleração possa ser considerada constante neste intervalo; todavia, intervalos muito pequenos fariam com que os tempos de cálculo da trajetória completa se tornassem proibitivos. Na prática, o Δt usado é da ordem de 0,5 a 1 fs (1 fs =10ˉ¹5 s).

O tempo da trajetória avaliada é em geral curto, visto que a simulação de dinâmica é um processo computacionalmente bastante custoso. Para sistemas na faixa dos milhares de átomos, a maior parte das simulações se restringe às centenas de picossegundos (1 ps = 10ˉ¹² s), a menos que computadores com elevada capacidade de processamento sejam utilizados. Para sistemas menores, a trajetória pode em alguns casos alcançar a faixa dos nanossegundos.

Há programas de domíno público para cálculos de dinâmica molecular como, por exemplo, o Tinker e o GROMACS.

2.2.3 Métodos ab initio

Os programas quanto‐mecânicos ab initio não incluem aproximações além das discutidas até aqui e representam os orbitais através dos chamados conjuntos de base.10,11 Um conjunto de base mínimo contém apenas um número de funções necessário para acomodar todos os elétrons de um átomo, por exemplo, uma função 1s para os átomos de hidrogênio e hélio, funções 1s e 2s para o lítio e o berílio e assim por diante. O maior problema do conjunto de base mínimo é a impossibilide de se expandir ou contrair os orbitais para que se ajustem ao ambiente molecular.

Uma maior flexibilidade na descrição dos elétrons é conseguida com os chamados conjuntos de base de valência dividida, nos quais as funções que representam os elétrons de valência são divididas em dois componentes, um mais interno e compacto e outro mais externo e difuso. A introdução de funções de número quântico secundário maior(funções do tipo p para átomos de hidrogênio e funções do tipo d para os demais) nos conjuntos de base de polarização permite que pequenos deslocamentos do centro de carga eletrônica em relação às posições nucleares sejam possíveis. Os conjuntos de base contendo funções difusas representam melhor sistemas aniônicos e estados excitados. Estes conjuntos são obtidos pela adição de orbitaiss e p muito difusos para melhorar a descrição de pares de elétrons de alta energia.

Há programas gratuitos que permitem cálculos ab initio com essas e outras funções de base, como o GAMESS.12

2.2.4 Deformação no comprimento de ligação

Se uma determinada ligação é comprimida ou estirada a energia sobe (Figura 4). Desta forma, o campo de força é parametrizado com valores de distância, por exemplo, para uma ligação C-C com hibridização sp3 e odesvio destes valores acarreta um aumento na energia da molécula.

Figura 4: Deformação no comprimento de ligação

A função potencial para a deformação angular (figura 5) deve considerar os diferentes tipos de átomos e hibridização molecular. Assim, o campo de força deve ser capaz de contemplar casos especiais, como por exemplo, o ciclobutano.

Figura 5: Deformação no ângulo de ligação

2.2.5 Barreira de energia de rotações intramoleculares (ângulos de torsão)

A análise conformacional envolve a rotação do ângulo de torsão θ que é formado por quatro átomos (A1, A2, A3 e A4). Estas rotações intramoleculares exigem energia. Na Figura 6, o valor de ângulo de torsão θ é de 180° e posiciona os átomos A1 e A4 o mais distante possível (conformação mais estável). Mudanças no ângulo de torsão θ ocasionam uma aproximação dos referidos átomos provocando um aumento de energia do sistema.

Figura 6:Ângulo de torsão (θ) com um valor de 180º

2.2.6 Interações de van der Waals

O raio de van der Waals de um átomo é o seu tamanho efetivo. Quando dois átomos não ligados são aproximados, a atração de van der Waals entre eles aumenta (decréscimo na energia). Quando a distância entre eles é igual à soma dos raios de van der Waals, a atração é máxima. Se os átomos são aproximados ainda mais, ocorre uma forte repulsão de van der Waals (Figura7).

Figura 7: Interações de van der Waals.

2.2.7 Métodos semi-empíricos

A equação de Schrödinger é uma equação diferencial e sua resolução envolve a avaliação de um grande número de integrais. No caso de cálculos com métodos ab initio, o número de integrais cresce aproximadamente com a quarta potência do número de funções de base, chegando a alguns milhões mesmo para moléculas pequenas.

Nos métodos quânticos semi‐empíricos, a negligência de um grande número dessas integrais foi a solução adotada para economizar tempo de máquina e também reduzir a quantidade de memória necessária nos cálculos. Um grande esforço foi dedicado na busca de métodos para a introdução de parâmetros empíricos ou previamente calculados na resolução da equação de Schrödinger que permitissem a eliminação do cálculo de algumas dessas integrais.

Após os trabalhos pioneiros de Pople e colaboradores, Dewar e colaboradores desenvolveram uma série de programas com o objetivo de tornar acessíveis cálculos semi‐empíricos de orbitais moleculares aos não especialistas. Além disso, os programas deveriam fornecer informações estruturais quimicamente precisas com um custo razoável de tempo de cálculo. Sem dúvida, os métodos semi‐empíricos de maior uso até hoje foram o AM1 (Austin Model 1 ou Modelo Austin 1)14 e o PM3 (Parametric Method 3 ou Método Paramétrico 3)15. Ambos os métodos incorporam aproximações semelhantes, mas diferem nas suas parametrizações.

Os métodos semi‐empíricos têm aplicação limitada a sistemas que contêm elementos para os quais foram desenvolvidos parâmetros correspondentes. No caso do método AM1, estes elementos são H, Na, K, Rb, Zn, Hg, B, Al, C, Si, Ge, Sn, N, P, O, S, F, Cl, Br e I; o método PM3 tem parâmetros para os mesmos elementos(exceto o B) e para outros elementos que não estão contidos no AM1: Be, Mg, Cd, Ga, In, Tl, As, Sb, Bi, Se e Te. A parametrização para estes métodosfoi desenvolvida para reproduzir uma série de dados experimentais, incluindo geometrias de equilíbrio, calores de formação, momentos de dipolo e energias de ionização. Alguns destes métodos estão presentes em pacotes de programas de domínio público, como o popular programa de cálculo de orbitais moleculares MOPAC,16 e em muitos pacotes comerciais. Recentemente, foram apresentadas duas novas parametrizações do modelo semi‐empírico, o RM1 (Recife Model 1)17 e o PM6 (Parametric Method 6)18, que disponibiliza parâmetros para 70 elementos da tabela periódica, incluindometais de transição.

2.2.8 Funcional de Densidade

Um modelo alternativo ao baseado em orbitais moleculares é o modelo do funcional de densidade (DFT). Nestes modelos, considera‐se que a energia de um conjunto de elétrons sob influência de um campo externo é um funcional único da densidade eletrônica. Esta dependência aparece em dois termos da energia eletrônica, chamados funcional de troca e funcional de correlação.19 Como a energia é expressa como uma função de uma única “variável”, a densidade eletrônica (que é função das 3 coordenadas cartesianas), as equações que resultam da aplicação deste modelo são mais simples do que as resultantes da teoria de Hartree‐Fock, onde as “variáveis” são o conjunto de funções de onda de um elétron (que são funções de 3N variáveis, onde N é o número de átomos do sistema).

Alguns funcionais foram desenvolvidos a partir da mecânica quântica fundamental e outros a partir da parametrização de funções que melhor reproduzem resultados experimentais. Desse modo, pode‐se dizer que há versões ab initio e semi‐empíricas do modelo DFT. Atualmente, um dos modelos mais utilizados é o modelo do funcional de troca híbrido de 3 parâmetros de Becke e do funcional de correlação de Lee‐Yang‐Parr (B3LYP), devido à qualidade dos seus resultados, particularmente para moléculas orgânicas.

2.2.9 Métodos comparativos

Além dos métodos de modelagem baseados na avaliação da energia, métodos de modelagem por comparação com estruturas conhecidas têm aplicações importantes para o estudo de compostos bioativos. É o caso da modelagem de proteínas por homologia, visto que em muitos casos somente a seqüência primária de um receptor ou de uma enzima é conhecida. A modelagem do receptor/enzima permite que estudos de interação com os bioligantes sejam aplicados da mesma forma que no caso de estruturas tridimensionais ‘reais’ de proteínas.

Estes métodos utilizam um algoritmo de reconhecimento de moldes adequados em um banco de dados de estruturas tridimensionais de proteínas, que depende do reconhecimento da similaridade entre as estruturas. Então, a partir de um ou mais moldes, os resíduos de aminoácidos da proteína em estudo são montados para se criar um modelo primário da estrutura, que é posteriormente otimizado por um método de mecânica molecular, em geral. Naturalmente, quanto mais resíduos forem idênticos entre o molde e o modelo, melhor será a qualidade do modelo final.

Este procedimento pode ser feito em sítios de domínio público disponíveis na Internet, como o servidor Swiss‐Model,21,22 ou através de programas específicos.

3 APLICAÇÃO DE MODELAGEM MOLECULAR EM QUÍMICA MEDICINAL

Neste tópico descrevem-se alguns resultados recentes de projetos desenvolvido por um grupo de pesquisas na aplicação de técnicas de modelagem molecular em Química Medicinal.

3.1 Estudos de Inibidores de Cicloxigenase (Cox) e 5-Lipoxigenase (5-LO): Proposta de novos inibidores Seletivos de 5- LO

O ácido araquidônico (ácido 5-Z,8-Z,11-Z,14-Z-icosatetraenóico, AA), liberado a partir de fosfolipídios de membranas celulares pela enzima fosfolipase A235,36, desempenha um papel importante na inflamação. Como substrato para as enzimas 5-Lipoxigenase (5-LO) e Cilcoxigenase (COx) - também conhecida como prostaglandina endoperóxido sinta se (PGHS) - o ácido araquidônico produz, respectivamente, o hidroperóxido (AA-OOH) que é subseqüentemente reduzido a leucotrienos (LTs), e o PGH2 que origina as prostaglandinas (PGs) (Fig. 8). Revisões recentemente publicadas descrevem vários aspectos dos inibidores de 5-LO37-39.

Figura 8 : Esquema do metabolismo de ácido araquidônico (AA)

O LTB4 é um potente agente quimiotáctil para as células inflamatórias (polimorfo nucleares, PMN) e pode ter um papel importante na última fase da inflamação, comumente observada em pacientes asmáticos. Os LTC4, LTD4 e LTE4, coletivamente identificados como substâncias de reação lenta da anafilaxia (SRS-A), têm efeitos farmacológicos potentes sobre as contrações dos músculos lisos, o estimulo da secreção de muco bronquial e o aumento da permeabilidade vascular, sendo considerados mediadores importantes dos processos alérgicos40,41.

O emprego de técnicas de modelagem molecular permitiu que fossem calculadas superfícies de energia potencial (parciais) utilizando o Hamiltoniano AM142-45 para vários derivados de pirazolinas (1-2) e de indazolinonas (3-7), descritos como inibidores de COx e 5-LO. Foram também estudados alguns derivados de 5-tioaril-N-fenil-pirazóis (8a-g) sintetizados no LASSBio46 que, embora estruturalmente relacionados com compostos anti-edematogênicos da classe dos 5-arilamino-N-fenil-pirazóis, não apresentaram atividade antiinflamatória significativa.. Estes estudos visaram a identificação de propriedades estereoeletrônicas, dependentes das estruturas, que pudessem estar relacionadas à atividade inibitória. A partir dos resultados e considerando o mecanismo molecular de ação correntemente aceito para este tipo de atividade farmacológica, foi possível propor novos inibidores mais seletivos da enzima 5-LO42-45.

A 5-LO é uma enzima ferro-não-heme dependente, encontrada primariamente em PMN e eosinófilos. Como dito acima, esta enzima catalisa a bioformação dos leucotrienos a partir do ácido araquidônico (Figura 8), estando envolvida, portanto, nos eventos inflamatórios40,47. O mecanismo de ação de Lipoxigenase proposto por Musser & Kreft48, Chasteen et. al.49 e Schilstra et. al.50 foi adotado em nossos estudos. Estes autores consideram a participação de um ciclo redox Fe(II)/Fe(III) neste mecanismo oxidativo.

Cucurou et. al.47 indicaram que no caso da inibição da Lipoxigenase-1 de soja (L-1), os inibidores BW755c (1) e fenidona (2) são oxidados a um cátion-radical que seria a espécie ativa na inibição da atividade enzimática. A etapa determinante do processo parece ser a oxidação enzimática do átomo de nitrogênio heterocíclico (N1), destes inibidores (1-7), pelo Fe(III)48. A Figura 9 ilustra o mecanismo proposto. O processo de inativação de L-1 inicia-se com a oxidação da fenidona a um cátion-radical (F+•)51. Este cátion-radical (F+•) ou outras espécies radicalares (F•) geradas pela atividade peroxidase de L-1 parecem estar envolvidos na inativação irreversível de Lipoxigenase, entre elas a L-1 (Figura 9).

Figura 9: Possível mecanismo molecular da inativação da Lipoxigenase de soja (L-1) por 7 e 8 proposto por Cucurou et. al.47

Ensaios com agentes redutores tais como o ácido ascórbico e derivados de tióis (glutationa=GSH; etanotiol=RSH) foram realizados47, objetivando analisar o papel das espécies radicalares no processo de inativação de L-1. Os resultados demonstraram que ocorre redução dos derivados radicalares (F• ou F+•) à forma neutra (F) (Fig. 9). Entretanto, a ação nucleofílica do ácido ascórbico e dos tióis não foi observada com os metabólitos eletrofílicos da fenidona já que estudos com fenidona [14C] não revelaram a presença de altos marcados47.

A oxidação posterior de F•• ou F+•a 1-fenil-1H-pirazol-3-ol poderia ocorrer através de diferentes reações (Fig.9). Primeiramente, por dismutação, gerando fenidona e seu derivado desidrogenado (Fig. 9, via a) 42. Em segundo lugar, a oxidação pela L-1 poderia levar a desidro-fenidona via carbocátion ou íon imônio, F+ (Fig. 9, via b). Esta espécie poderia formar ligação covalente com a proteína por reação com um resíduo nucleofílico de um aminoácido. Esta hipótese tem suporte experimental na correlação linear observada entre a formação de ligação covalente e a inativação da lipoxigenase47. O papel do oxigênio na inativação da L-1 sugere, finalmente, a possibilidade de oxidação do cátion-radical (F+•) a desidro-fenidona e O2•- (via c). Esta última espécie (O2•-) poderia ser responsável pela inativação da L-1 por um mecanismo relacionado à oxidação da metionina (Fig. 9).

Na Tabela I estão descritos os dados físico-químicos mais relevantes, calculados por AM1, assim como as atividades farmacológicas destes derivados. Os compostos 1-7 (Fig. 10, Tabela I) são inibidores52 de COx e 5-LO, sendo o derivado ICI207968 (7)52 o mais seletivo a nível de 5-LO. Os compostos 8a-g (Fig. 10, Tabela I) foram sintetizados no LASSBio46 como prováveis bioisósteros com atividade antiinflamatória, apresentando, entretanto, fraca atividade. Os dados da Tabela I não evidenciaram qualquer relação entre os momentos de dipolo ou os potenciais de ionização e a atividade inibitória. Entretanto, as densidades eletrônicas (DE) sobre N1 são muito diferentes para os compostos ativos e inativos. Esta diferença sugere que altas densidades eletrônicas são desejáveis para a atividade.

Figura 10. Inibidores de 5-LO e COx (1-7 e 8ª-g)

As densidades eletrônicas de 8a-g parecem insuficientes para participarem do ciclo redox inibitório. Estes resultados suportam a hipótese mecanística mencionada anteriormente e justificam a fraca atividade observada para os derivados sintéticos da série 8. Os dados da Tabela I sugerem, também, que as DEs sobre N2 também são muito baixas para participarem do ciclo redox. A análise dos dados mostra que existe inversão de atividade através da 5-LO e COx para os compostos 2 (DE N1=0,54) e 3 (DE N1=0,51), o que parece sugerir que valores menores de DE (3-7) estão associados com maior seletividade para 5-LO.

Como uma extensão desta hipótese, uma série de derivados triazólicos (9a-f) com as características estruturais mínimas necessárias à atividade desejada foram planejados, considerando o perfil farmacológico procurado42. Estes compostos, derivados do di-hidro-1,2,3-triazol (9a-f) (Fig. 11, Tabela II), são análogos estruturais do BW755c (1)47 e da fenidona (2) e possuem um segundo anel fenila de forma a potencializar suas propriedades hidrofóbicas e, em consequência, aumentar a seletividade frente a 5-LO53.

Figura 11. Inibidores proposto 9a-f para a 5-LO (derivados 1,2,3-triazólicos)

O padrão de substituição em C5 foi concebido de forma a permitir estudos mecanísticos. A II descreve os resultados obtidos para estes compostos pelo método AM1. Considerando a hipótese descrita acima, que correlaciona a atividade com a densidade eletrônica, espera-se que os derivados 9a, 9e e 9f se comportem como o composto 7, inibindo seletivamente a enzima 5-LO. Espera-se, por outro lado, que os compostos 9c e 9d sejam inativos.

Foi proposta, ainda, uma segunda série de derivados heterocíclicos baseada nos mesmos critérios farmacológicos, os derivados 10a-h (Fig. 12, Tabela III) 45. A Tabela III contém resultados obtidos pelo método AM1. Como os compostos 8a-g (Tabela I), as DEs sobre N1 destes derivados são muito baixas, impedindo a participação no ciclo redox Fe(II)/Fe(III). Entretanto, algumas densidades eletrônicas sobre N2 tem os valores desejados, o que permite antecipar que o composto 10e deve agir como um inibidor de ambas as enzimas 5-LO e COx, sendo seletivo frente a 5-LO. O composto 10a, no qual se observou uma baixa densidade eletrônica sobre N2, deve ser ativo em relação a 5-LO, mas inativo em relação à enzima COx.

Figura 12. Inibidores propostos 10a-h para 5- LO (derivados pirrólicos e pirazólicos)

Os resultados obtidos nestes estudos47 antecipam a hipótese de que a modulação da carga residual sobre o átomo de nitrogênio heteroaromático pode determinar as propriedades inibidoras sobre COx/5-LO, permitindo o planejamento de inibidores seletivos de 5-LO. Outrossim, a conjugação de propriedades hidrofóbicas adequadas, determinará o perfil farmacocinético que as novas substâncias poderão apresentar, viabilizando sua avaliação farmacológica in vivo.

3.2 Estudo Teórico da Tromboxana A2: Proposta de uma Conformação Bioativa

Como mencionado anteriormente, a Cilcoxigenase (COx), também conhecida como prostaglandina endoperóxido sintase (PGHS), é responsável pela ciclização-oxidativa do ácido araquidônico (AA, 11) a endoperóxido de prostaglandina H2 (PGH2, 12), que é isomerizado pela Tromboxana sintase (TXS) à Tromboxana A2 (TXA2, 13), uma substância com meia vida de 3 minutos nas condições fisiológicas e que se transforma em Tromboxana B2 (TXB2, 14), um metabólito estável com perfil biológico35 distinto (Fig. 13).

Figura 13. Biotransformação da TXA2 e da TXB2 (14) a partir do ácido araquidônico (11)

A TXA2 (13) é um vasoconstritor potente no sistema cardiovascular e um indutor potente da agregação plaquetária e da reação de liberação plaquetária no sistema sangüíneo, atuando a nível de receptores específicos. O PGH2 (12) possui um perfil farmacológico semelhante ao da TXA2 (13) com propriedades agonísticas no receptor de TXA235. Nestes sistemas biológicos a TXB2 (14) é inativa. O receptor de TXA2(TP) é comumente referido como receptor de TXA2 e PGH2 uma vez que estes eicosanóides produzem efeitos biológicos semelhantes. De fato, análogos estáveis do endoperóxido natural PGH2 (12) tais como U-44609 (15) e U-46619 (16) (Fig. 14) mimetizam os efeitos da TXA254 (13).

Figura 14. Estrutura das prostanóides U-46609 (15) e U-46619 (16). Análogos estáveis de PGH2 (12) e TXA2 (13)

Recentemente, a seqüência de aminoácidos do receptor TP humano foi deduzida a partir da seqüência de nucleotídeos do clone de cDNA que o codifica, utilizando-se tecido de placenta humana e de um clone parcial, obtido a partir de uma cultura de células leucêmicas megacariocíticas humanas55. Acredita-se que o receptor TP pertença à família de receptores acoplados à proteína-G (GPCR) por apresentar sete domínios transmembrânicos, uma característica dos GPCR, e uma significativa homologia seqüencial com o receptor de rodopsina, um membro desta família55,56.

Estes conhecimentos permitiram que o receptor de TXA2 fosse modelado56 por mecânica molecular, com base em sua seqüência de aminoácidos e na estrutura tridimensional da bacteriorodopsina, obtida por criomicroscopia eletrônica. Esta enzima também contém sete segmentos transmembrânicos, mas como apresenta pouca homologia seqüencial com a rodopsina56, o modelo gerado deve ser usado com cautela, na falta da estrutura cristalográfica.

O sítio do receptor assim modelado apresenta dois resíduos de aminoácidos, Arg295 e Ser201, separados por uma região hidrofóbica. A Arg295 parece interagir com a carboxila terminal da TXA255 (13) De fato, o ajuste da TXA2 (13) ao receptor modelado permite identificar as interações entre o grupamento ácido carboxílico terminal da TXA2 (13) e o resíduo de Arg295, e entre a hidroxila alílica em C-15 e o resíduo de Ser201.

Estes estudos, porém, não indicaram nenhuma interação envolvendo os átomos de oxigênio do sistema bis-oxabiciclo[3.1.1]-heptano e outros resíduos de aminoácido do receptor (Fig. 15).

Figura 15. Representação esquemática do sítio receptor da Tromboxana A2(TXA2, 13) a partir do receptor TP modelado por Yamamoto et. al56, demonstrando a interação entre o receptor TP e a TXA2 (13)

A estrutura da TXB2 (14) (Fig. 13), o metabólito estável da TXA2 (13) que não possui propriedades agregantes plaquetária, determinada por difração de raios-X57, mostrou duas formas de cristalização, denominadas α e β (Fig. 16, a e b) de conformação semelhante, exceto quanto à orientação das cadeias α58, com a função ácido carboxílico terminal orientada em direções opostas. Neste estudo, foi sugerido que a estrutura da TXB2 (14) pode ser comparada com a da TXA2 (13) porque, apesar da TXB2 (14) ser inativa a nível das plaquetas, a preferência de ligação dos anéis endocíclicos são suficientemente diferentes para assegurar um reconhecimento molecular distinto entre estes dois prostanóides57.

Figura 16. Estruturas cristalográficas α (a) e β (b) da TXB2 (14)57

O estudo de diversos agonistas e antagonistas do receptor de TXA2 (13) por FTIR em solução diluída de CCl4, mostrou a formação de ligação hidrogênio intramolecular em muitos deles59,60. Por exemplo, o composto U-46619 (16) (Fig. 14), um conhecido eficiente agonista do receptor TP, apresenta 80% das suas moléculas na conformação estabilizada por ligação hidrogênio intramolecular entre a carboxila e a hidroxila, formando um macrociclo de quinze membros, relativamente rígido59. Nestes estudos, o CCl4 foi usado como solvente porque sua constante dielétrica (ε=2,23 a 25oC) é semelhante à encontrada no interior das proteínas. Os autores concluíram que embora a conformação predominante não deva ser a bioativa, porque compromete a função carboxila, considerada um farmacoforo importante para as atividades agonista e antagonista, esta não deve diferir muito daquela determinada em CCl459,60. Este tipo de argumento não considera, entretanto, que o microambiente em que o agonista ou o antagonista interagem no receptor, isto é, o sítio receptor, é formado por resíduos de aminoácidos diversos, que podem inclusive atuar como doadores ou aceptores em ligações hidrogênio, alterando a conformação da molécula. Este microambiente, portanto, não é perfeitamente mimetizado por um solvente apolar incapaz de interagir com a molécula por ligações hidrogênio.

Pelo contrário, é razoável supor que compostos contendo, simultaneamente, grupos doadores e aceptores de hidrogênio, diluídos em solventes apolares, formem, preferencialmente, ligações hidrogênio intramoleculares, desde que a geometria molecular o permita, e o fator entrópico não seja desfavorável. De qualquer forma, como salientado no próprio trabalho59,60, o estudo conformacional da TXA2 (13) em solução diluída de CCl4 pode ser útil como subsídio para a modelagem da conformação através de cálculos teóricos, usualmente realizados com a molécula isolada, simulando o vácuo (ε=1) 59,60.

Ezumi et. al.61 também estudaram por mecânica molecular e orbitais moleculares (MNDO) a TXA2 (13) e o composto U-46619 (16) dentre outros compostos, e assumiram que a conformação bioativa não deve diferir muito das conformações mais estáveis no estado gasoso. Estes autores propuseram duas conformações bioativas para estes dois prostanóides, uma envolvendo ligação hidrogênio intramolecular, semelhante à encontrada em CCl4, e outra, mantendo basicamente a geometria anterior, porém, com a carboxila do ácido terminal e a hidroxila em C-15, opostamente orientadas61.

Recentemente, foi proposto por mecânica molecular e orbitais moleculares um modelo espacial do sítio farmacofórico do receptor de Tromboxana A262, com base no estudo das conformações mais estáveis de cinco conhecidos antagonistas de receptor de TXA2.

Nossos resultados sobre a análise conformacional da TXA2 (13) pelo método AM144,45,63,64, conseguiram classificar as conformações obtidas para a TXA2 (13) em dois grupos distintos. Um, representado pela conformação 1 (Fig. 17 a), onde ocorre formação de ligação hidrogênio intramolecular, como no caso do análogo U-46619 (16)59,61, e o outro, representado pela conformação 2 (Fig. 17, b), na qual as cadeias α e ω58 estão relativamente afastadas e, tanto a carboxila do ácido terminal como a hidroxila em C-15 estão orientadas em direções opostas, porém distintas daquelas propostas por Ezumi et. al.61.

Figura 17. Confirmações I (a) e 2 (b) da TXA2 obtidas por AMI45

Estes resultados foram obtidos a partir de um estudo de sobreposição molecular, no qual foram considerados o átomo de carbono C-1 da carboxila terminal, um dos átomos de oxigênio do anel endoperóxido ligado em C-9, e o átomo de oxigênio da hidroxila em C-15 da estrutura da TXA2 (13) (Fig. 13), para a sobreposição dos pares de conformações em análise.

As conformações 1 e 2 da TXA2 (13), selecionadas por representarem conformações distintas, foram sobrepostas, então, às conformações α e β da TXB2 (14) obtidas por difração de raios-X57. As sobreposições da conformação α da TXB2 (14) com as conformações 1 e 2 da TXA2 (13) apresentaram desvios médios de 3,3 e 0,8 Å, respectivamente, enquanto que as sobreposições da conformação β da TXB2(14) com as conformações 1 e 2 da TXA2 (13) apresentaram, desvios médios de 2,7 e 0,6 Å, respectivamente (Tabela IV).

Assim, as sobreposições obtidas com os melhores níveis de similaridade conformacional, ocorreram entre as conformações α ou β da TXB2 (14) e a conformação 2 da TXA2 (13). Cabe ressaltar que a TXB2 (14) não apresenta atividade em nível de receptor de TXA2 (13) e, portanto, a melhor sobreposição poderia ser desfavorável para a conformação 2 da TXA2 (13). Entretanto, observando-se melhor a estrutura da TXB2, pode-se notar que o sistema oxacíclico apresenta duas hidroxilas, conferindo elevada hidrofilia a esta região da molécula, distintamente ao que ocorre na TXA2 (13). Como conseqüência, este autacóide não deve ser reconhecido pelo receptor TP, muito provavelmente, devido a esta diferença de caráter lipofílico ao nível da subunidade estrutural complementar ao sítio hidrofóbico estereo-exigente do receptor TP62 e não porque apresentaria uma conformação inadequada.

Os resultados destes estudos indicam que a conformação 2 pode representar a conformação bioativa da TXA2 (13), possuindo como referência, uma conformação similar da TXB2 (14), no estado sólido, podendo ser então utilizada como modelo farmacofórico no planejamento de novos antagonistas de receptor TP45,63.

4 MÉTODOS DE MODELAGEM PARA ESTUDO E PLANEJAMENTO DE COMPOSTOS BIOATIVOS

A aplicação de métodos computacionais no estudo e no planejamento de compostos bioativos tem se tornado uma prática rotineira nos dias atuais. Do ponto de vista do planejamento, é importante destacar que, quando se fala do uso da modelagem molecular, não se pretende chegar a uma molécula bioativa simplesmente através do uso de programas de computador. O processo de desenvolvimento dessas moléculas, devido à sua complexidade, envolve necessariamente o trabalho de uma equipe multidisciplinar, que emprega um grande conjunto de métodos computacionais de modo sistemático de forma a facilitar e otimizar o processo de desenvolvimento de compostos bioativos, em uma constante troca de informações com grupos de síntese química e avaliação da atividade destes compostos.

Esses métodos computacionais podem ser usados como ferramentas do planejamento racional de compostos bioativos, assim chamado porque é orientado por uma hipótese racional sobre o mecanismo de ação destes compostos.1 A ação das moléculas bioativas é um fenômeno bastante complexo, mas um dos paradigmas da química medicinal é que essas moléculas têm seus efeitos associados a interações ou até reações químicas com estruturas macromoleculares presentes no sistema vivo, proteínas, na sua grande maioria.1 Caso estas proteínas sejam receptores celulares, as moléculas bioativassão classificadas como agonistas ou antagonistas de receptores; no caso das enzimas, estas moléculas atuam como inibidores enzimáticos.

5 CONCLUSÃO

O número de programas que permitem o uso da modelagem molecular no estudo e no planejamento de compostos bioativos continua aumentando todos os anos. Contudo, quase todos dependem da determinação da estrutura e da energia moleculares e, portanto,têm um ou mais dos métodos discutidos nesse trabalho como ponto de partida. Como escolher o melhor método de modelagem para um determinado problema? De modo simples, a escolha de um ou de outro método depende da

consideração de três fatores principais: as quantidades que se deseja determinar, a precisão desejada (necessária) e a capacidade de cálculo disponível. Por exemplo, um estudo conformacional simples pode normalmente ser feito com um modelo de mecânica molecular, que é bastante rápido e apresenta resultados de boa qualidade, mas quando se está interessado em quantidades eletrônicas, como a densidade eletrônica ou a energia dos orbitais, ou então no estudo de reações químicas, onde ligações são quebradas e formadas, a escolha obrigatoriamente passa por um método quântico. Mas a precisão com que se quer que essas quantidades sejam obtidas determina qual o método quântico de escolha. Métodos semi‐empíricos, em geral, produzem resultados menos precisos do que os métodos ab initio, por exemplo.

Contudo, o uso de métodos mais sofisticados, necessariamente computacionalmente mais custosos, pode serinviável dependendo do tamanho do sistema molecular e das características do equipamento disponível, como a memória e a velocidade do processador.

REFERÊNCIAS BIBLIOGRÁFICAS

Abraham, D. J. X-ray crystallografy and drug design. In: Perun, T. J. & Propst, C. L. (eds) Computer-aided drug design: methods and applications. Marcel Dekker, New York, 1989

Academic Press, 1996. 361p.

SANT’ANNA, C.M.R. Glossário de termos usados no planejamento de fármacos (recomendações da IUPAC para 1997). Quim. Nova, v.25, n.3, p.505-512, 2002.

PARTINGTON et al. Molecular Pharmacology, v.8, p. 269-270, 1972.

HÖLTJE, H. D.; FOLKERS, G. Molecular Modeling- Basic principles and applications. NewYork: VHC, 1996. 194 p

Burkert, U. & Allinger, N. L. Molecular mechanics. ACS Monograph 177. American Chemical Society, Washington, DC, 1982

Charretier, E. & Guéron, M. Biochem. Cell. Biol. 1991, 69, 322.

Campbell, W. B. Lipid-derived autacoids: eicosanoids and platelet-activating factor. In: Gilman, A. G., Rall, T. W., Nies, A. S. & Taylor, P. (eds) Goodman & Gilman's the pharmacological basis of therapeutics. 8th ed. Pergamon Press, New York, 1990.

Campbell, S. F. Química Nova 1991, 14, 196.

Barreiro, E. J. e Fraga, C. A. M.; Química Medicinal, As Bases Moleculares da Ação dos Fármacos, Artmed Editora, Porto Alegre, 2001

Leach, A. R. Molecular Modelling, Principles and Applications. Longman, Edinburgh, 1998.

Stewart, J. J. P. MOPAC: A Semiempirical MolecularOrbital Program. ESCOM, Leiden, 1990.