【量子AIが変える医療 #5】第2章 QAOA創薬最適化 ハミルトニアン構築と最適化実行

第2章：量子近似最適化アルゴリズム（QAOA）を用いた創薬最適化

QAOAは組み合わせ最適化問題を量子コンピューターで解くアルゴリズムです。タンパク質-リガンド結合の最大化問題をQUBO（Quadratic Unconstrained Binary Optimization）に変換して解きます。

from qiskit.circuit.library import QAOAAnsatz
from qiskit.quantum_info import SparsePauliOp
from qiskit.primitives import Estimator
from scipy.optimize import minimize

class DrugMoleculeOptimizer:
    """
    QAOAを用いた薬物分子の結合部位最適化
    タンパク質-リガンド結合エネルギー最大化問題をQUBOに変換
    """

    def __init__(self, n_sites: int, p: int = 2):
        """
        Args:
            n_sites: 結合候補部位数
            p: QAOAの深さパラメータ
        """
        self.n_sites = n_sites
        self.p = p
        self.estimator = Estimator()

    def build_interaction_hamiltonian(
        self,
        binding_energies: np.ndarray,
        interaction_matrix: np.ndarray
    ) -> SparsePauliOp:
        """
        結合エネルギーと相互作用からハミルトニアンを構築

        Args:
            binding_energies: 各部位の結合エネルギー [n_sites]
            interaction_matrix: 部位間相互作用行列 [n_sites, n_sites]
        Returns:
            コスト関数ハミルトニアン
        """
        n = self.n_sites
        pauli_list = []

        # 線形項（結合エネルギー）
        for i, energy in enumerate(binding_energies):
            pauli_str = 'I' * (n - 1 - i) + 'Z' + 'I' * i
            pauli_list.append((pauli_str, -energy / 2))

        # 2次項（部位間相互作用）
        for i in range(n):
            for j in range(i + 1, n):
                if abs(interaction_matrix[i, j]) > 1e-6:
                    pauli_str = ['I'] * n
                    pauli_str[n - 1 - i] = 'Z'
                    pauli_str[n - 1 - j] = 'Z'
                    pauli_list.append(
                        (''.join(pauli_str), interaction_matrix[i, j] / 4)
                    )

        return SparsePauliOp.from_list(pauli_list)

    def optimize(
        self,
        binding_energies: np.ndarray,
        interaction_matrix: np.ndarray,
        n_iterations: int = 100
    ) -> dict:
        """
        QAOA最適化の実行

        Returns:
            最適解の情報（選択部位、推定結合エネルギー）
        """
        hamiltonian = self.build_interaction_hamiltonian(
            binding_energies, interaction_matrix
        )
        qaoa_circuit = QAOAAnsatz(
            cost_operator=hamiltonian,
            reps=self.p
        )
        n_params = len(qaoa_circuit.parameters)
        initial_params = np.random.uniform(0, 2 * np.pi, n_params)

        def cost_function(params):
            job = self.estimator.run(
                [qaoa_circuit], [hamiltonian], [params]
            )
            return job.result().values[0]

        result = minimize(
            cost_function,
            initial_params,
            method='COBYLA',
            options={'maxiter': n_iterations}
        )

        return {
            'optimal_energy': result.fun,
            'optimization_success': result.success,
            'n_iterations': result.nit
        }