トポロジカル極性表面積
学生時分には、球や円錐の「表面積」の計算を教わりましたが、
何処で使う計算なのか把握していなかったこともあって、
個人的に、随分とぞんざいに扱ってきた気がしています。
特に、ほぼ同時に習う「体積」と比較して、
高校範囲の物理・化学では幽霊のような存在だったことが大きな要因です。
円周2πr --（積分）--> 円の表面積4πr^2 --（積分）--> 円の体積(3/4)πr^3
を直感的に理解できれば見通しが俄然よくなるのですが、
大学物理で場の理論を押さえていないと、やはり味気ないのです。
すみません、話がそれました。
極性表面積（Polar Surface Area, PSA）とは、その名のとおり、
分子の表面のうち、極性を帯びている部分の面積です。
分子表面の極性が高いと、物理的に細胞膜を通過しにくい、
すなわち、吸収されにくいと考えられます。
さて、計算ですが、
面積なので、それぞれの分子の三次元構造を前提に積分計算するのが
まっとうな印象ですが、欠点として、計算に時間がかかる点、そして
正確な三次元座標の生成が要求される点が挙げられます。
（まともな三次元構造を生成するのは、けっこう難しいものです）
「トポロジカル極性表面積（TPSA）」は、SMILES形式のファイルから
手軽に極性表面積を近似計算する方法です。
その考え方はXLogPと似ていて、回帰分析により、
43種類の原子タイプに係数を割り振るというものです。
Ertlらの原著論文*1に具体的な値が表にまとめられています。
CDKでは、以下のように実装されています。
/* tpsa.java */import java.io.*;
import org.openscience.cdk.io.iterator.IteratingSMILESReader;
import org.openscience.cdk.interfaces.IMolecule;
import org.openscience.cdk.qsar.IMolecularDescriptor;import org.openscience.cdk.CDKConstants;
import org.openscience.cdk.Ring;
import org.openscience.cdk.annotations.TestClass;
import org.openscience.cdk.annotations.TestMethod;
import org.openscience.cdk.aromaticity.CDKHueckelAromaticityDetector;
import org.openscience.cdk.exception.CDKException;
import org.openscience.cdk.interfaces.IAtom;
import org.openscience.cdk.interfaces.IAtomContainer;
import org.openscience.cdk.interfaces.IBond;
import org.openscience.cdk.interfaces.IRingSet;
import org.openscience.cdk.qsar.DescriptorSpecification;
import org.openscience.cdk.qsar.DescriptorValue;
import org.openscience.cdk.qsar.IMolecularDescriptor;
import org.openscience.cdk.qsar.result.DoubleResult;
import org.openscience.cdk.qsar.result.IDescriptorResult;
import org.openscience.cdk.ringsearch.AllRingsFinder;
import org.openscience.cdk.tools.CDKHydrogenAdder;
import org.openscience.cdk.tools.manipulator.AtomContainerManipulator;
import java.util.ArrayList;
import java.util.HashMap;
import java.util.List;class tpsa {
public static void main(String args[]){
boolean checkAromaticity = false;
HashMap map = new HashMap();// positions in atom profile strings are: symbol, max-bond-order, bond-order-sum,
// number-of-neighbours, Hcount, formal charge, aromatic-bonds, is-in-3-membered-ring,
// single-bonds, double-bonds, triple-bonds.map.put("N+1.0+3.0+3+0+0+0+0+3+0+0", new Double(3.24)); // 1
map.put("N+2.0+3.0+2+0+0+0+0+1+1+0", new Double(12.36)); // 2
map.put("N+3.0+3.0+1+0+0+0+0+0+0+1", new Double(23.79)); // 3
map.put("N+2.0+5.0+3+0+0+0+0+1+2+0", new Double(11.68)); // 4
map.put("N+3.0+5.0+2+0+0+0+0+0+1+1", new Double(13.6)); // 5
map.put("N+1.0+3.0+3+0+0+0+1+3+0+0", new Double(3.01)); // 6
map.put("N+1.0+3.0+3+1+0+0+0+3+0+0", new Double(12.03)); // 7
map.put("N+1.0+3.0+3+1+0+0+1+3+0+0", new Double(21.94)); // 8
map.put("N+2.0+3.0+2+1+0+0+0+1+1+0", new Double(23.85)); //9
map.put("N+1.0+3.0+3+2+0+0+0+3+0+0", new Double(26.02)); // 10
map.put("N+1.0+4.0+4+0+1+0+0+4+0+0", new Double(0.0)); //11
map.put("N+2.0+4.0+3+0+1+0+0+2+1+0", new Double(3.01)); //12
map.put("N+3.0+4.0+2+0+1+0+0+1+0+1", new Double(4.36)); //13
map.put("N+1.0+4.0+4+1+1+0+0+4+0+0", new Double(4.44)); //14
map.put("N+2.0+4.0+3+1+1+0+0+2+1+0", new Double(13.97)); //15
map.put("N+1.0+4.0+4+2+1+0+0+4+0+0", new Double(16.61)); //16
map.put("N+2.0+4.0+3+2+1+0+0+2+1+0", new Double(25.59)); //17
map.put("N+1.0+4.0+4+3+1+0+0+4+0+0", new Double(27.64)); //18
map.put("N+1.5+3.0+2+0+0+2+0+0+0+0", new Double(12.89)); //19
map.put("N+1.5+4.5+3+0+0+3+0+0+0+0", new Double(4.41)); //20
map.put("N+1.5+4.0+3+0+0+2+0+1+0+0", new Double(4.93)); //21
map.put("N+2.0+5.0+3+0+0+2+0+0+1+0", new Double(8.39)); //22
map.put("N+1.5+4.0+3+1+0+2+0+1+0+0", new Double(15.79)); //23
map.put("N+1.5+4.5+3+0+1+3+0+0+0+0", new Double(4.1)); //24
map.put("N+1.5+4.0+3+0+1+2+0+1+0+0", new Double(3.88)); //25
map.put("N+1.5+4.0+3+1+1+2+0+1+0+0", new Double(14.14)); //26map.put("O+1.0+2.0+2+0+0+0+0+2+0+0", new Double(9.23)); //27
map.put("O+1.0+2.0+2+0+0+0+1+2+0+0", new Double(12.53)); //28
map.put("O+2.0+2.0+1+0+0+0+0+0+1+0", new Double(17.07)); //29
map.put("O+1.0+1.0+1+0+-1+0+0+1+0+0", new Double(23.06)); //30
map.put("O+1.0+2.0+2+1+0+0+0+2+0+0", new Double(20.23)); //31
map.put("O+1.5+3.0+2+0+0+2+0+0+0+0", new Double(13.14)); //32map.put("S+1.0+2.0+2+0+0+0+0+2+0+0", new Double(25.3)); //33
map.put("S+2.0+2.0+1+0+0+0+0+0+1+0", new Double(32.09)); //34
map.put("S+2.0+4.0+3+0+0+0+0+2+1+0", new Double(19.21)); //35
map.put("S+2.0+6.0+4+0+0+0+0+2+2+0", new Double(8.38)); //36
map.put("S+1.0+2.0+2+1+0+0+0+2+0+0", new Double(38.8)); //37
map.put("S+1.5+3.0+2+0+0+2+0+0+0+0", new Double(28.24)); //38
map.put("S+2.0+5.0+3+0+0+2+0+0+1+0", new Double(21.7)); //39map.put("P+1.0+3.0+3+0+0+0+0+3+0+0", new Double(13.59)); //40
map.put("P+2.0+3.0+3+0+0+0+0+1+1+0", new Double(34.14)); //41
map.put("P+2.0+5.0+4+0+0+0+0+3+1+0", new Double(9.81)); //42
map.put("P+2.0+5.0+4+1+0+0+0+3+1+0", new Double(23.47)); //43if(args.length != 1){
System.err.println("tpsa <smiles-file>");
System.exit(1);
}
FileInputStream fis = null;
IteratingSMILESReader isr = null;
try{
fis = new FileInputStream(new File(args[0]));
isr = new IteratingSMILESReader(fis);
} catch (Exception e) {
e.printStackTrace();
}while( isr.hasNext() ){ // Read a file line by line
IMolecule imol = (IMolecule)isr.next();
String id = (String)imol.getProperty("cdk:Title"); // The 2nd columnIRingSet rs;
List<String> profiles = new ArrayList<String>();
IAtomContainer ac;
ac = (IAtomContainer)imol;
double tpsa = 0.0;try{
rs = (new AllRingsFinder()).findAllRings(ac);
if (checkAromaticity) {
AtomContainerManipulator.percieveAtomTypesAndConfigureAtoms(ac);
CDKHueckelAromaticityDetector.detectAromaticity(ac);
}// iterate over all atoms of ac
java.util.Iterator atoms = ac.atoms().iterator();
while ( atoms.hasNext() ) {
IAtom atom = (IAtom) atoms.next();
if ( ! atom.getSymbol().equals("N") &&
! atom.getSymbol().equals("O") &&
! atom.getSymbol().equals("S") &&
! atom.getSymbol().equals("P") ) { continue; }
int singleBondCount = 0;
int doubleBondCount = 0;
int tripleBondCount = 0;
int aromaticBondCount = 0;
double maxBondOrder = 0;
double bondOrderSum = 0;
int hCount = 0;
int isIn3MemberRing = 0;// counting the number of single/double/triple/aromatic bonds
List<IBond> connectedBonds = ac.getConnectedBondsList(atom);
for (IBond connBond : connectedBonds) {
if (connBond.getFlag(CDKConstants.ISAROMATIC))
aromaticBondCount++;
else if (connBond.getOrder() == CDKConstants.BONDORDER_SINGLE)
singleBondCount++;
else if (connBond.getOrder() == CDKConstants.BONDORDER_DOUBLE)
doubleBondCount++;
else if (connBond.getOrder() == CDKConstants.BONDORDER_TRIPLE)
tripleBondCount++;
}
int formalCharge = atom.getFormalCharge();
java.util.List connectedAtoms = ac.getConnectedAtomsList(atom);
int numberOfNeighbours = connectedAtoms.size();// EXPLICIT hydrogens
for (int nei = 0; nei < numberOfNeighbours; nei++)
if ( ( (IAtom) connectedAtoms.get(nei) ).getSymbol().equals("H") )
hCount++;
// IMPLICIT hydrogens
Integer implicitHAtoms = atom.getHydrogenCount();
if (implicitHAtoms == CDKConstants.UNSET) {
implicitHAtoms = 0;
}for (int hi = 0; hi < implicitHAtoms; hi++) {
hCount++;
numberOfNeighbours++;
singleBondCount++;
}
// Calculate bond order sum using the counters of 1/2/3 bonds
bondOrderSum += singleBondCount * 1.0;
bondOrderSum += doubleBondCount * 2.0;
bondOrderSum += tripleBondCount * 3.0;
bondOrderSum += aromaticBondCount * 1.5;
// setting maxBondOrder
if (singleBondCount > 0)
maxBondOrder = 1.0;
if (aromaticBondCount > 0)
maxBondOrder = 1.5;
if (doubleBondCount > 0)
maxBondOrder = 2.0;
if (tripleBondCount > 0)
maxBondOrder = 3.0;// isIn3MemberRing checker
if (rs.contains(atom)) {
IRingSet rsAtom = rs.getRings(atom);
for (int ri = 0; ri < rsAtom.getAtomContainerCount(); ri++) {
Ring ring = (Ring) rsAtom.getAtomContainer(ri);
if (ring.getRingSize() == 3)
isIn3MemberRing = 1;
}
}
// create a profile of the current atom (atoms[atomIndex])
String profile = atom.getSymbol() + "+" + maxBondOrder +
"+" + bondOrderSum + "+" +
numberOfNeighbours + "+" + hCount + "+" + formalCharge +
"+" + aromaticBondCount + "+" +
isIn3MemberRing + "+" + singleBondCount + "+" +
doubleBondCount + "+" + tripleBondCount;
//logger.debug("tpsa profile: "+ profile);
profiles.add(profile);
}
// calculate the tpsa for the AtomContainer ac
for (int pi = 0; pi < profiles.size(); pi++) {
if (map.containsKey(profiles.get(pi))) {
tpsa += (Double) map.get(profiles.get(pi));
}
}
profiles.clear();
} catch (Exception e1) {
e1.printStackTrace();
}
System.out.printf("%s\t%.3f\n", id, tpsa);
}
}
}
このように、表面積を直接計算するのが難しいときは、
他の性質から近似的に求めていくことが実用上可能です。
大事なことは、表面積を求めることに意味があるか、じっくり吟味することだと思います。
