Cleiton Breder Eller

• Pós-Graduação em Ecologia, Universidade Estadual de Campinas

• ellercb@gmail.com

Condições abióticas formam um importante filtro que pode determinar a ocorrência de organismos em uma dada região (Lambers et al 2008). Os limites de tolerância a condições abióticas que organismos apresentam são em grande parte determinados por suas respostas fisiológicas a fatores como temperatura, humidade, radiação solar, etc.

Plantas apresentam uma grande diversidade de respostas fisiológicas a condições ambientais adversas. Para ilustrar essa diversidade podemos usar a diversidade de respostas fisiológicas que plantas podem apresentar durante períodos de seca. Algumas plantas perdem total ou parcialmente suas folhas, diminuindo gastos metabólicos durante o período de seca e retomam suas atividades fisiológicas normais no final do período adverso. Outras plantas podem manter suas folhas e evitar esse período de “hiato” fisiológico, mas para preservar a integridade de seus tecidos essas plantas adotam estratégias mais conservadoras em seu uso de recursos, como por exemplo, diminuindo a perda de água foliar durante certos períodos do dia através do fechamento estomático. Entretanto, algumas plantas podem apresentar tecidos mais resistentes a escassez de água do que outras e dessa forma apresentar uma forma menos restritiva de controle estomático durante períodos secos.

As relações entre respostas fisiológicas de plantas a condições ambientais frequentemente assumem formas não-lineares, especialmente quando plantas são expostas a condições abióticas severas e processos fisiológicos operam próximo de um ponto de colapso. Um exemplo disso é a função sigmoidal ou exponencial que descreve a perda de condutividade do xilema de plantas quando submetidas a estresse hídrico (Cochard et al 2013), ou a função, frequentemente aceita, como hiperbólica que descreve a regulação estomática de folhas em resposta a seca atmosférica (Damour et al 2010).

Logo, uma abordagem possível para explicar e predizer respostas fisiológicas de plantas a condições abióticas é através do ajuste de modelos aos dados usando o método de máxima verossimilhança. Muitas vezes diversos modelos alternativos podem ser ajustados ao mesmo conjunto de dados. Um método adequado para selecionar o modelo que melhor se ajusta aos dados é através da comparação de modelos alternativos usando o Critério de Informação de Akaike (AIC).

O ajuste de modelos usando o método de máxima verossimilhança baseia-se em encontrar o parâmetro que maximiza a função de verossimilhança mantendo constante o conjunto de dados observados. A função de verossimilhança para uma única observação pode ser escrita como:

Onde é a função de densidade probabilística que indica a verossimilhança de que , que representa uma variável aleatória, seja igual a , o valor observado. O valor de máxima verossimilhança indica o valor de parâmetro ( ) que tornam o valor observado mais verossímil. Para se calcular a verossimilhança de um conjunto de dados basta-se calcular o produto das funções de verossimilhança de cada observação independente (Batista, 2009).

Embora o valor de máxima verossimilhança possa ser encontrado analiticamente para modelos simples, algoritmos de otimização numérica são empregados para estimar o valor de máxima verossimilhança na superfície de verossimilhança gerada por modelos complexos (Purcell, 2007).

O uso da razão de verossimilhança na comparação de modelos com um número diferente de parâmetros pode ser um método enviesado, pois modelos com mais parâmetros tendem a apresentar valores maiores de verossimilhança do que modelos com menos parâmetros. Esse enviesamento viola o princípio de parcimônia, o que pode gerar modelos com baixo poder de generalização para novos conjuntos de dados.

O critério de informação de Akaike (AIC) é frequentemente usado na comparação de modelos com números de parâmetros diferentes, pois penaliza o número de parâmetros presentes nos modelos a serem comparados:

Onde é a verossimilhança máxima do modelo e é o número de parâmetros do modelo. Como o AIC é baseado na distância relativa de Kullback-Leibler entre um modelo e o modelo verdadeiro (que é desconhecido), maiores valores de AIC indicam uma maior distância relativa entre o modelo em questão e o modelo verdadeiro. O modelo com menor valor de AIC dentre um conjunto de modelos é aquele que mais se aproxima do modelo verdadeiro que gerou os dados. Como o AIC é uma medida de distância relativa, ele n seu valor absoluto não possui significado. Logo, mesmo o modelo com menor valor de AIC pode ser uma representação ruim dos dados, se o conjunto de modelos iniciais for composto somente por modelos que não se adequam aos dados (Burnham & Arnderson, 2002).

Nesse exemplo, eu usarei o método de máxima verossimilhança para ajustar 5 modelos alternativos que explicam a relação entre a condutância estomática da copa (gs; mmol H2O m-2 s-1) de Drimys brasiliensis (Miers) e déficit de pressão de vapor atmosférico (DPV). Em seguida eu selecionarei o modelo mais plausível usando o AIC.

Os dados de gs usados nesse exemplo foram estimados a partir de dados de fluxo de seiva (medidos com o método de razão de calor; Burgess et al 2001) do tronco de cinco árvores, usando uma equação de Penman-Montheth modificada (Oren et al 2001). Dessa forma, esses dados integram o comportamento estomático de toda a copa da árvore. Os dados de déficit de pressão de vapor atmosférico (DPV; kPa) foram calculados usando a equação descrita em Allen et al. (1998) e representam a umidade necessária para saturar completamente o ar a uma dada temperatura, de forma que valores mais altos indicam uma atmosférica mais seca.

data<-data.frame(gs= c(141.176646, 194.331243, 98.163473, 76.590556, 82.529139, 87.606701, 59.452335, 42.558477,      40.115762,  38.311903, 38.388689,  37.666063,  35.906439,  41.474234,  15.463164,  18.194682,  20.337496,    21.406662,  22.957422,  19.833383, 20.815754,  20.455995,  18.589086,  18.208495,  18.348146,  17.462362,  16.160288,  15.187332,  16.727786,  16.360037, 22.092856,  20.376686,  19.399405,  22.012389,  19.708153,  22.586931,  20.750567,  18.760256,  17.148380,  26.884228, 27.396807,  31.305529,  33.307537,  34.793092,  34.338288,  33.664394,  18.912528,  13.209450,  16.689383,  18.513572, 16.753321,  15.562446,  18.316990,  16.440499,  15.692933,  14.871380,  19.755290,  29.372969,  39.208676, 30.919732,  32.579054,  44.356648,  14.221484,   8.468941,  12.007191,  12.138791,  10.441959,  10.749526,   9.996734, 15.289351,  21.494993,  13.390422,  10.678711,  12.425899,  14.290464,  14.874564,  18.402725,  19.915527), DPV = c(0.05581582,0.05680164, 0.17358928, 0.26819740, 0.27948776, 0.28406637, 0.44355298, 0.60238018, 0.69014739, 0.73810399,0.70979785,0.66042682, 0.61454498, 0.45251073, 0.78253424, 0.87402125, 1.00969468, 1.08642385, 1.10307149, 1.37480871, 1.40175124, 1.46593474, 1.64865285, 1.68875295, 1.65933791, 1.71204056, 1.83471854, 1.78161821, 1.56631407, 1.45427959, 0.30104750, 0.56446840, 0.64639925, 0.83745318, 1.03393778, 1.03784102, 1.14086899, 1.26688868, 1.37043206, 0.87721769, 0.84388881, 0.64798411, 0.56082325, 0.48349315, 0.44798315, 0.41469726, 0.57365967, 1.04014364, 1.02292313, 1.05722667, 1.18283286, 1.39152029, 1.17751143, 1.28852221, 1.17828942, 1.27627788, 0.84631404, 0.61313355, 0.42253148, 0.45310876, 0.34935628, 0.25734419, 0.68070389, 1.41664728, 1.21342971, 1.32789002, 1.64890874, 1.67011601, 1.77204579, 1.03450311, 0.81865940, 1.16508390, 1.13995946, 1.08710844, 1.04327017, 0.94798441, 0.70208438, 0.66802293) )

O DPV é a força motriz do fluxo transpiratório em plantas. Como altos fluxos transpiratórios podem causar danos ao sistema de transporte de água de plantas (cavitação) plantas tendem a limitar o fluxo transpiratório através da diminuição de gs com o aumento do DPV (Damour et al 2010).

Monteith (1995) propõe uma relação linear (1) entre gs e DPV, descrita pela função:

gs = gmax – s*DPV (1)

Onde o intercepto da função (gmax) é o maior valor de gs alcançado pela planta e a inclinação (s) representa sensibilidade estomática ao aumento de DPV. Embora esse modelo preveja valores de gs negativos em situações de alto DPV, Monteth sugere que tais valores dificilmente são encontrados em situações naturais. Entretanto, a função hiperbólica (2) proposta anteriormente por Lohammar et al (1980) evita esse problema ao fixar um limiar mínimo para gs com o parâmetro D:

gs = gmax/(1+DPV/D) (2)

Franks (2014) sugere que a maior parte dos trabalhos relacionando gs com DPV produz funções próximas de um decaimento exponencial (3) ou uma função logarítmica negativa (4), logo também testarei esses dois tipos de modelos:

gs = gmax*sDPV (3)

gs= -s*ln(DPV) (4)

Todos estes modelos predizem uma resposta imediata do gs a incrementos de DPV. A fim de testar a existência de um período inicial de tolerância do gs a incrementos de DPV usarei uma função de decaimento sigmoidal (5) onde o parâmetro g50 representa o ponto em que o gs alcançará 50% do valor inicial:

gs = gmax/(1+g50*sDPV) (5)

Eu construi funções para o cálculo da log-verossimilhança negativa de cada modelo. Em seguida estimei os parâmetros que minimizavam a log-verossimilhança negativa dessas funções usando a função mle do pacote bbmle:

#modelo linear (1)
LL.m1<-function(gmax,s,sigma){
  m=gmax-s*data$DPV
  -sum(dnorm(data$gs,mean=m,sd=sigma,log=T))
}
m1<-mle(LL.m1,start=list(gmax=100,s=50,sigma=50))
summary(m1)
#modelo hiperbólico (2)
LL.m2<-function(gmax,D,sigma){
  m=gmax/(1+(data$DPV/D))
  -sum(dnorm(data$gs,mean=m,sd=sigma,log=T))
}
m2<-mle(LL.m2,start=list(gmax=200,D=10,sigma=50))
summary(m2)
#modelo de decaimento exponencial (3)
LL.m3<-function(gmax,s,sigma){
  m=gmax*s^data$DPV
  -sum(dnorm(data$gs,mean=m,sd=sigma,log=T))
}
m3<-mle(LL.m3,start=list(gmax=200,s=10,sigma=50))
summary(m3)
#modelo log negativo (4)
LL.m4<-function(s,sigma){
  m=-s*log(data$DPV)
  -sum(dnorm(data$gs,mean=m,sd=sigma,log=T))
}
m4<-mle(LL.m4,start=list(s=10,sigma=50))
summary(m4)
#modelo de decaimento sigmoidal (4)
LL.m5<-function(gmax,g50,s,sigma){
  m=gmax/(1+g50*s^data$DPV)
  -sum(dnorm(data$gs,mean=m,sd=sigma,log=T))
}
m5<-mle(LL.m5,start=list(gmax=200,g50=50,s=5,sigma=50))
summary(m5)

A comparação entre os modelos foi feita usando a função AICtab do pacote bbmle:

AICtab(m1,m2,m3,m4,m5,weights=T)
     dAIC  df weight
m2   0.0 3  1     
m3  56.4 3  <0.001
m5  73.9 4  <0.001
m1 121.5 3  <0.001
m4 136.3 2  <0.001

E por fim, podemos comparar graficamente a performance dos modelos em relação aos dados observados:

plot(gs~DPV,data,ylim=c(-25,200),xlab="VPD(kPa)",ylab="gs(mmol m-2 s-1)")
abline(0,0)
curve(coef(m1)[1]-coef(m1)[2]*x, from=0, to=2, add=T)
curve(coef(m2)[1]/(1+(x/coef(m2)[2])), from=0, to=2, add=T,col="red")
curve(coef(m3)[1]*coef(m3)[2]^x, from=0, to=2, add=T,col="blue")
curve(-coef(m4)[1]*log(x), from=0, to=2, add=T,col="green")
curve(coef(m5)[1]/(1+coef(m5)[2]*coef(m5)[3]^x),from=0, to=2, add=T,col="orange")
legend("topright", c("linear","hiperbólico","decaimento exponencial", "log negativo", "decaimento sigmoidal"),
       lty=c(1,1,1,1,1),col=c("black","red","blue","green","orange"))

Entre os modelos testados, o modelo que melhor se ajusta aos dados é a função hiperbólica de Lohammer et al (1980). A afirmação de que plantas não ocorrem em ambientes em que o DPV seja alto suficiente para fazer com que o modelo linear de Monteith (1995) faça predições de valores de gs negativos não foi corroborada em nosso exemplo. O modelo linear previa valores de gs negativo quando o DPV fosse maior que ~1.7 kPa e em nossas observações o DPV chegou até 1.83 kPa. Além disso, o modelo de Monteith não fez predições precisas em valores baixos de DPV (menores que 0.5 kPa). Drimys brasiliensis também não demonstrou nenhuma evidência de um período de tolerância inicial ao aumento de DPV em nosso exemplo (dado a rejeição do modelo de decaimento sigmoidal). Um próximo passo para essa aplicação seria incluir outras variáveis ambientais no modelo como disponibilidade água no solo e radiação fotossinteticamente ativa.

• Allen, R. G., Pruitt, W. O., Wright, J. L., Howell, T. A., Ventura, F., Snyder, R., … & Elliott, R. (2006). A recommendation on standardized surface resistance for hourly calculation of reference ETo by the FAO56 Penman-Monteith method. Agricultural Water Management, 81(1), 1-22.
• Batista, J.L.F(2009) Verossimilhança e Máxima Verossimilhança.
• Burgess, S. S., Adams, M. A., Turner, N. C., Beverly, C. R., Ong, C. K., Khan, A. A., & Bleby, T. M. (2001). An improved heat pulse method to measure low and reverse rates of sap flow in woody plants. Tree Physiology, 21(9), 589-598.
• Burnham K.P., Anderson D.R. (2002). Model Selection and Multimodel Inference A Practical Information-Theoretic Approach. Springer-Verlag, New York
• Cochard, H., Badel, E., Herbette, S., Delzon, S., Choat, B., & Jansen, S. (2013). Methods for measuring plant vulnerability to cavitation: a critical review. Journal of experimental botany, ert193.
• DAMOUR, G., SIMONNEAU, T., COCHARD, H. and URBAN, L. (2010), An overview of models of stomatal conductance at the leaf level. Plant, Cell & Environment, 33: 1419–1438. doi: 10.1111/j.1365-3040.2010.02181.x
• Franks, P. J. (2004). Stomatal control and hydraulic conductance, with special reference to tall trees. Tree Physiology, 24(8), 865-878.
• Lambers, H., Chapin, F., Pons, T.L. (2008) Plant Physiological Ecology, , Springer-Verlag, New York
• Lohammer T., Larsson S., Linder S. & Falk S.O. (1980) FAST_ Simulation models of gaseous exchange in Scots pine. Ecological Bulletin 32, 505–523.
• Monteith J.L. (1995) A reinterpretation of stomatal responses to humidity. Plant, Cell & Environment 18, 357–364.
• Oren, R., Sperry, J. S., Ewers, B. E., Pataki, D. E., Phillips, N., & Megonigal, J. P. (2001). Sensitivity of mean canopy stomatal conductance to vapor pressure deficit in a flooded Taxodium distichum L. forest: hydraulic and non-hydraulic effects. Oecologia, 126(1), 21-29.
• Purcell S. (2007) Maximum likelihood primer. http://statgen.iop.kcl.ac.uk/bgim/mle/sslike_1.html

Este ensaio é um produto de disciplina da pós-graduação da Universidade de São Paulo. Para citá-lo:

Eller, CB. 2014. Modelando respostas fisiológicas de plantas em função de condições microclimáticas. In: Prado , P.I & Batista, J.L.F. Modelagem Estatística para Ecologia e Recursos Naturais. Universidade de São Paulo. url: http://cmq.esalq.usp.br/BIE5781.