预留度数Rx表达式的推导

1. sympy简介

  sympy是Python中强大的符号运算包，可用于多种计算，如微积分运算、解方程、矩阵运算等

  人工晶状体度数的计算代码参考R程序包enbrown/iol-calculations ¹

  安装该包的方法

install.packages("devtools")
library(devtools)
install_github("enbrown/iol-calculations")

  该R程序包包含了公开的多种计算公式，提供了每种公式中ELP的计算，继而从ELP推导出人工晶状体的度数。但是作者并没有提供确定了P之后预留度数Rx的计算代码

  以下利用sympy演示Holladay1及HofferQ计算公式下Rx表达式的推导

2. Holladay1

2.1 P的计算源码

enbrown/iol-calculations包内源代码：

#' Holladay 1 Formula for Effective Lens Position
#' 
#' Calculate IOL effective lens position for emmetropia given axial length, corneal radius of curvature, 
#' and lens constant.
#' 
#' Note: If the Holladay Surgeon-factor constant for the lens is not provided, it is calculated from
#' the lens A-constant or pACD-constant. Similarily, if the corneal radius of curvature is not provided,
#' it is calculated from the corneal curvature (in diopters) and corneal index of refraction.
#' 
#' @param L length of the eye in millimeters (mm)
#' @param R corneal radius of curvature (mm)
#' @param K corneal power (D) used to determine corneal radius of curvature
#' @param cornea_n corneal index of refraction used to determine corneal radius of curvature from K
#' @param S IOL surgeon-factor constant (mm)
#' @param A IOL A constant (D) used to determine surgeon-factor
#' @param pACD IOL personalized ACD constant (mm) used to determine surgeon-factor
#' @return ELP in mm
#' @seealso \code{\link{ELP}}
#' @family ELP
#' @references \url{https://encrypted.google.com/books?id=NhWJsGFK6qgC&pg=PA8}
Holladay.1.ELP <- function(L, R, K, cornea_n, S, A, pACD) {
  args <- list(L = L)
  #https://encrypted.google.com/books?id=NhWJsGFK6qgC&pg=PA8&lpg=PA8&dq=colenbrander+formula&source=bl&ots=j__sC0XHHg&sig=DEaK9Qg0LOvedzq4itcfe7pBZ4Q&hl=en&sa=X&ei=scdlUuC1JcHP2wW90IDAAg&ved=0CD4Q6AEwAg#v=onepage&q=colenbrander%20formula&f=false
  if (missing(R) || ! is.finite(R)) {
    if (missing(cornea_n) || !is.finite(cornea_n)) {
      cornea_n <- Constants$Biometry$corneal_index
      warning("Using standard corneal index of refraction to convert K to R for Holladay 1 equation: ", cornea_n)
    }
    args$cornea_n <- cornea_n
    args$K <- K
    R <- 1000 * (cornea_n - 1) / K
    warning("Convert K to R for Holladay 1 equation: ", R, "mm")
  } else {
    args$R <- R
  }
  if (missing(S) || ! is.finite(S)) {
    if (! missing(A) && is.finite(A)) {
      S <- Constants$IOL$A_to_S_a0 + Constants$IOL$A_to_S_a1 * A
      args$A <- A
      warning("Using Holladay approximation for S from A constant for Holladay 1 equation: ", S, "mm")
    } else if (! missing(pACD) && is.finite(pACD)) {
      S <- Constants$IOL$pACD_to_S_a0 + Constants$IOL$pACD_to_S_a1 * pACD
      args$pACD <- pACD
      warning("Using Holladay approximation for S from pACD constant for Holladay 1 equation: ", S, "mm")
    } else {
      stop("No S (surgeon factor) provided for Holladay 1 equation")
    }
  } else {
    args$S <- S
  }
  AG <- L * 12.5 * 1 / 23.45
  temp <- R * R - AG * AG / 4
  if (temp < 0) {
    warning("Calculation of anatomic anterior chamber depth poorly defined for this combination of corneal curvature and axial length in Holladay 1 formula")
    temp <- 0
    result <- NA
    attr(result, 'parameters') <- args
    return(result)
  }
  aACD <- 0.56 + R  - sqrt(temp)
  ELP <- aACD + S
  result <- ELP
  attr(result, 'parameters') <- args
  return(result)
}
ELP.functions$Holladay.1 <- Holladay.1.ELP

#' Holladay-1 Formula for IOL Power
#' 
#' Calculate IOL power for emmetropia given axial length, 
#' corenal curvature, and effective lens position.
#' 
#' Note: if the corneal radius, \code{R}, is not provided, it is determined 
#' from the corneal power, \code{K}, and the corneal index of refraction, 
#' \code{cornea_n}.
#' 
#' @param L length of the eye in millimeters (mm)
#' @param R corneal radius of curvature (mm)
#' @param ELP IOL effective lens position (mm)
#' @param K corneal power (D) used to determine corneal radius of curvature
#' @param cornea_n corneal index of refraction used to determine corneal radius of curvature from K, defaults to 1.3375
#' @param aqueous_n aqueous index of refraction, defaults to 1.3375
#' @param RT retinal thickness (mm), defaults to 0.2 mm
#' @param Rx desired resulting refractive error (D), defaults to 0 D
#' @param V vertex distance of Rx (mm), defaults to 13 mm
#' @return Power of IOL (D)
#' @seealso \code{\link{Power}}
#' @family Power
#' @references \url{https://encrypted.google.com/books?id=NhWJsGFK6qgC&pg=PA8}
Holladay.1.Power <- function(L, R, ELP,
                                       K, cornea_n = 1.3375,
                                       aqueous_n =1.336, RT = 0.2, Rx = 0, V = 13) {
  args <- list(L = L)
  #https://encrypted.google.com/books?id=NhWJsGFK6qgC&pg=PA8&lpg=PA8&dq=colenbrander+formula&source=bl&ots=j__sC0XHHg&sig=DEaK9Qg0LOvedzq4itcfe7pBZ4Q&hl=en&sa=X&ei=scdlUuC1JcHP2wW90IDAAg&ved=0CD4Q6AEwAg#v=onepage&q=colenbrander%20formula&f=false
  if (missing(R) || ! is.finite(R)) {
    if (missing(cornea_n) || !is.finite(cornea_n)) {
      cornea_n <- Constants$Biometry$corneal_index
      warning("Using standard corneal index of refraction to convert K to R for Holladay 1 equation: ", cornea_n)
    }
    args$cornea_n <- cornea_n
    args$K <- K
    R <- 1000 * (cornea_n - 1) / K
    warning("Convert K to R for Holladay 1 equation: ", R, "mm")
  } else {
    args$R <- R
  }
  Alm <- L + RT
  nc <- 4/3
  naRnc1Alm <- aqueous_n * R - (nc - 1) * Alm
  naRnc1ELP <- aqueous_n * R - (nc - 1) * ELP
  P <- (1000 * aqueous_n * (naRnc1Alm - 0.001 * Rx * (V * naRnc1Alm + Alm * R))) /
    ((Alm - ELP) * (naRnc1ELP - 0.001 * Rx * (V * naRnc1ELP + ELP * R)))
  result <- P
  attr(result, 'parameters') <- args
  return(result)
}
Power.functions$Holladay.1 <- Holladay.1.Power

  计算P所需的参数：

L, R, ELP,K, cornea\_n = 1.3375,aqueous\_n =1.336, RT = 0.2, Rx = 0, V = 13

  K转化为R：
$$R=\frac{1000\times (cornea\_n-1)}{K}$$

  调整眼轴等：
$$Alm=L+RT$$
$$nc=\frac{4}{3}$$

$$naRnc1Alm=aqueous\_n\times R-(nc-1)\times Alm$$
$$naRnc1ELP=aqueous\_n\times R-(nc-1)\times ELP$$
$$AG=\frac{12.5\times L}{23.45}$$

  ELP的表达式：
$$ELP=0.56+R-\sqrt{(}{R}^2-\frac{A{G}^2}{4})+S$$

  将ELP计算结果代入下式即可计算人工晶体度数P：

$$P=\frac{1000\times aqueous\_n\times (naRnc1Alm-0.001\times Rx\times (V\times naRnc1Alm+Alm\times R))}{(Alm-ELP)\times (naRnc1ELP-0.001\times Rx\times (V\times naRnc1ELP+ELP\times R))}$$

2.2 Rx表达式的推导

  由上述表达式从P反求Rx，运算小心一点可能也可以算出来Rx，哈哈~~
  不过终究还是复杂了些。这时候我们借助sympy将会十分省力，使用python console，我们可以交互式输入

#####Holliday1
###ELP
from sympy import *         ##导入sympy
K,L,S = symbols('K,L,S')    ##申明一些基本符号
R = 1000 * (1.3375 - 1) / K
AG = L * 12.5 * 1 / 23.45
temp = R * R - AG * AG / 4
ELP = 0.56 + R  - sqrt(temp)+S
ELP

  于是我们得到ELP的表达式：

S - 337.5*sqrt(-6.23629694574426e-7*L**2 + K**(-2)) + 0.56 + 337.5/K

  即：
$$ELP=S-337.5\times \sqrt{(}-6.23629694574426e-7\times L^2+(1/K{)}^2)+0.56+337.5/K$$
  注：6.23629694574426e为科学计数法

上述输出的ELP表达式我们可以粘贴到记事本或者sublime之后，查找替换对应参数为Excel表格即可用于Excel内的计算，注意一下Excel内公式的表达与此代码的差异并修正即可

  有了ELP之后我们来导出Rx的表达式

from sympy import *       
P0,L,K,nc,Rx,ELP = symbols('P0,L,K,nc,Rx,ELP')
R = 1000 * (1.3375 - 1) / K
Alm = L + 0.2
naRnc1Alm = 1.336 * R - (nc - 1) * Alm
naRnc1ELP = 1.336 * R - (nc - 1) * ELP
P = (1000 * 1.336 * (naRnc1Alm - 0.001 * Rx * (13 * naRnc1Alm + Alm * R))) /((Alm - ELP) * (naRnc1ELP - 0.001 * Rx * (13 * naRnc1ELP + ELP * R)))
Rx =solve(P-P0, Rx)[0]  #解方程P-P0=0,Rx等于什么，此处用P0替代P
Rx

  Rx的表达式为：

1000.0*(-50.0*ELP**2*K*P0*nc + 50.0*ELP**2*K*P0 + 50.0*ELP*K*L*P0*nc - 50.0*ELP*K*L*P0 + 10.0*ELP*K*P0*nc - 10.0*ELP*K*P0 + 22545.0*ELP*P0 - 66800.0*K*L*nc + 66800.0*K*L - 13360.0*K*nc + 13360.0*K - 22545.0*L*P0 - 4509.0*P0 + 30120120.0)/(-650.0*ELP**2*K*P0*nc + 650.0*ELP**2*K*P0 + 16875.0*ELP**2*P0 + 650.0*ELP*K*L*P0*nc - 650.0*ELP*K*L*P0 + 130.0*ELP*K*P0*nc - 130.0*ELP*K*P0 - 16875.0*ELP*L*P0 + 289710.0*ELP*P0 - 868400.0*K*L*nc + 868400.0*K*L - 173680.0*K*nc + 173680.0*K - 293085.0*L*P0 + 22545000.0*L - 58617.0*P0 + 396070560.0)

  即：
$$Rx=1000.0\times \frac{-50.0\times EL{P}^2\times K\times P0\times \frac{4}{3}+50.0\times EL{P}^2\times K\times P0+50.0\times ELP\times K\times L\times P0\times \frac{4}{3}-50.0\times ELP\times K\times L\times P0+10.0\times ELP\times K\times P0\times \frac{4}{3}-10.0\times ELP\times K\times P0+22545.0\times ELP\times P0-66800.0\times K\times L\times \frac{4}{3}+66800.0\times K\times L-13360.0\times K\times \frac{4}{3}+13360.0\times K-22545.0\times L\times P0-4509.0\times P0+30120120.0}{-650.0\times EL{P}^2\times K\times P0\times \frac{4}{3}+650.0\times EL{P}^2\times K\times P0+16875.0\times EL{P}^2\times P0+650.0\times ELP\times K\times L\times P0\times \frac{4}{3}-650.0\times ELP\times K\times L\times P0+130.0\times ELP\times K\times P0\times \frac{4}{3}-130.0\times ELP\times K\times P0-16875.0\times ELP\times L\times P0+289710.0\times ELP\times P0-868400.0\times K\times L\times \frac{4}{3}+868400.0\times K\times L-173680.0\times K\times \frac{4}{3}+173680.0\times K-293085.0\times L\times P0+22545000.0\times L-58617.0\times P0+396070560.0}$$
  同样Rx的表达式也可以转换为Excel计算公式，结合上一步ELP的计算即可算出Rx

3 HofferQ

3.1 P的计算源码

enbrown/iol-calculations包内源代码：

#' Hoffer-Q Formula for Effective Lens Position
#' 
#' Calculate IOL effective lens position for emmetropia given axial length, 
#' corenal curvature, and lens A-constant or pACD constant.
#' 
#' Note: If the pACD constant, \code{pACD} for the lens is not provided, it is calculated from
#' the lens A-constant, \code{A}, using the Holladay approximation.
#' 
#' @param L length of the eye in millimeters (mm)
#' @param K average corneal curvature of the eye in diopters (D)
#' @param pACD IOL personalized ACD constant (mm)
#' @param A IOL A constant (D) used to determine pACD
#' @return ELP in mm
#' @seealso \code{\link{ELP}}
#' @family ELP
Hoffer.Q.ELP <- function(L, K=337.5/R, R, pACD, A) {
  args <- list(L = L, K = K)
  if (missing(pACD) || ! is.finite(pACD)) {
    pACD <- ELP.functions$Holladay(A)
    args$A <- A
    warning("Using Holladay approximation for Hoffer Q's pACD from A constant: ", pACD, "mm")
  } else {
    args$pACD <- pACD
  }
  if (L > 23) {
    M <- -1
    G <- 23.5
  } else {
    M <- +1
    G <- 28
  }
  if (L < 18.5) L <- 18.5
  if (L > 31) L <- 31
  ELP <- pACD +
    0.3 * (L - 23.5) +
    (tan(K * pi / 180))^2 +
    ((0.1 * M * (23.5 - L)^2) *
       tan(0.1 * (G - L)^2 * pi / 180)) -
    0.99166
  result <- ELP
  attr(result, 'parameters') <- args
  return(result)
}
ELP.functions$Hoffer.Q <- Hoffer.Q.ELP

#' Hoffer-Q Formula for Emmetropic IOL Power
#' 
#' Calculate IOL power for emmetropia given axial length, 
#' corenal curvature, and effective lens position.
#' 
#' @param L length of the eye in millimeters (mm)
#' @param K average corneal curvature of the eye in diopters (D)
#' @param ELP IOL effective lens position (mm)
#' @param Rx resulting refractive error (D), defaults to 0 D
#' @param V resulting refractive error vertex distance (mm), defaults to 13 mm
#' @return Power of emmetropic IOL (D)
#' @seealso \code{\link{Power}}
#' @family Power
Hoffer.Q.Power <- function(L, K=337.5/R, R, ELP, Rx = 0, V = 13) {
  args <- list(L = L, K = K, ELP = ELP, Rx = Rx, V = V)
  R <- Rx / (1 - 0.012 * Rx)
  P <- 1336 / (L - ELP - 0.05) -
    (1.336 / (
      (1.336 / (K + R)) - ((ELP + 0.05) / 1000)))
  attr(P, 'parameters') <- args
  return(P)
}
Power.functions$Hoffer.Q <- Hoffer.Q.Power

  计算P所需参数：

L, K=337.5/R, R, pACD, A

$$\{\begin{array}{ll}M=-1& \text{if }L>23\\ G=23.5& \end{array}\{\begin{array}{ll}M=+1& \text{else }\\ G=28& \end{array}$$

$$L=\{\begin{array}{ll}18.5& \text{if }L<18.5\\ 31& \text{if }L>31\end{array}$$

$$R=\frac{Rx}{1-0.012\times Rx}$$

  ELP的表达式直接有了，为：

ELP = pACD +0.3 * (L - 23.5) +(tan(K * pi / 180))^2 +((0.1 * M * (23.5 - L)^2) *tan(0.1 * (G - L)^2 * pi / 180)) -0.99166

  即：
$$ELP=pACD+0.3\times (L-23.5)+(tan\frac{K\pi }{180}{)}^2+0.1\times M\times (23.5-L{)}^2\times tan\frac{0.1\times (G-L{)}^2\times \pi }{180}-0.99166$$

  由此用下式计算P：
$$P=\frac{1336}{L-ELP-0.05}-\frac{1.336}{\frac{1.336}{K+R}-\frac{ELP+0.05}{1000}}$$

3.2 Rx表达式的推导

from sympy import *
P0,L,ELP,K,Rx = symbols('P0, L,ELP,K,Rx')
R = Rx / (1 - 0.012 * Rx)
P = 1336 / (L - ELP - 0.05) -(1.336 / ((1.336 / (K + R)) - ((ELP + 0.05) / 1000)))
Rx =solve(P-P0, Rx)[0]
Rx

  可得Rx为：

Rx= 250.0*(-400.0*ELP**2*K*P0 + 400.0*ELP*K*L*P0 - 40.0*ELP*K*P0 + 534400.0*ELP*P0 + 20.0*K*L*P0 - 534400.0*K*L - K*P0 - 534400.0*L*P0 + 26720.0*P0 + 713958400.0)/(-1200.0*ELP**2*K*P0 + 100000.0*ELP**2*P0 + 1200.0*ELP*K*L*P0 - 120.0*ELP*K*P0 - 100000.0*ELP*L*P0 + 1613200.0*ELP*P0 + 60.0*K*L*P0 - 1603200.0*K*L - 3.0*K*P0 - 1608200.0*L*P0 + 133600000.0*L + 80410.0*P0 + 2141875200.0)

  即：
$$Rx=250.0\times \frac{-400.0\times EL{P}^2\times K\times P+400.0\times ELP\times K\times L\times P-40.0\times ELP\times K\times P+534400.0\times ELP\times P+20.0\times K\times L\times P-534400.0\times K\times L-K\times P-534400.0\times L\times P+26720.0\times P+713958400.0}{-1200.0\times EL{P}^2\times K\times P+100000.0\times EL{P}^2\times P+1200.0\times ELP\times K\times L\times P-120.0\times ELP\times K\times P-100000.0\times ELP\times L\times P+1613200.0\times ELP\times P+60.0\times K\times L\times P-1603200.0\times K\times L-3.0\times K\times P-1608200.0\times L\times P+133600000.0\times L+80410.0\times P+2141875200.0}$$

---

1. enbrown/iol-calculations 源码链接: https://github.com/enbrown/iol-calculations ↩︎