XPC, as a DNA repair tool, contributes to the maintenance of genomic integrity. Given its critical function, it is imperative that the expression of XPC by our gene circuits be meticulously regulated to avoid underexpression, which may result in functional deficiencies, or overexpression, which could precipitate adverse consequences.

Based on our gene dynamics modeling results, a key determinant of the final XPC expression is the translation efficiency—the amount of protein produced per unit of mRNA. In synthetic biology, an effective strategy for regulating mRNA translation efficiency is sequence optimization. In our design, the XPC-P2A-RUP2 mRNA transcribed from the engineered plasmid contains only the coding sequence (CDS). To maintain the amino acid sequence of the expressed protein and preserve its biological function, synonymous codon replacement is needed within the CDS. This ensures that the final protein expression level remains within a desired range without altering the protein's structure or function.

However, the mechanism by which synonymous codon substitution affects translation efficiency is complex. To simplify, we can conceptualize that the amount of protein a given mRNA can ultimately produce before degradation depends on two primary factors: how fast it translates (translation rate) and how long it lasts (stability). Therefore, the total protein produced from a single mRNA molecule before degradation can be represented as the integral of the translation rate over the duration of its existence. Mathematically, this can be expressed as:

\[ P = \int_0^T r(t) \, dt \]

Where:

• \( P \): the total amount of protein produced by a single mRNA.
• \( r(t) \): the rate of translation at time \( t \).
• \( T \): the lifespan of the mRNA (i.e., how long it exists before degradation).

The translation rate \([r(t)]\) represents the ribosome scanning speed on mRNA. When the tRNA corresponding to a codon is more abundant, the ribosome typically scans that region more quickly. There exists a bias in the usage of synonymous codons, either within a single mRNA or across the entire transcriptome. Transcripts that preferentially use codons of higher usage frequency or corresponding to more abundant tRNAs generally exhibit faster translation rates^[7]. Consequently, several widely used metrics have been developed to assess this bias, such as the Codon Adaptation Index (CAI) and the tRNA Adaptation Index (tAI).

It seems plausible to modify metrics such as the Codon Adaptation Index (CAI) or the tRNA Adaptation Index (tAI) by replacing synonymous codons to adjust translation rates and thus regulate final protein expression levels. However, recent studies have demonstrated that translation rate \([r(t)]\) can influence mRNA stability \((T)\) in a codon-dependent manner^[9]. Additionally, synonymous codon substitutions alter the mRNA secondary structure, which simultaneously affects both translation rate \([r(t)]\) and stability \((T)\)^[1]. Generally, a more compact mRNA structure decreases ribosome scanning speed and prolongs mRNA half-life.

We have explored several commonly used mRNA design methods and encountered the following issues:

• Insufficient modeling of the intricate relationship between mRNA sequences and translation efficiency (as seen in CAI-based approaches)
• Difficulty in capturing long-range dependencies within mRNA sequences (e.g., with traditional machine learning models), particularly given that our XPC-P2A-RUP2 transcript is 3993 nucleotides (nt) in length
• The optimization objective usually being increasing translation efficiency (e.g., LinearDesign^[10]), whereas we may aim to downregulate protein expression
• Non-coding regions included, not dedicated to CDS design (e.g., RNA-FM^[2])

Here, we employ deep neural networks to achieve end-to-end prediction from mRNA CDS sequence to protein expression level, and propose a solution towards next-generation CDS sequence design.

Our model CodonBERTER (CodonBERT with Elongated Range), is adapted from the pretrained LLM for mRNAs, CodonBERT^[5]. CodonBERT uses codons as token embeddings instead of nucleotides which may learn better representations for protein related downstream tasks. CodonBERT was pretrained with over 10 million mRNA CDS sequences from a diverse set of organisms. A codon is composed of 3-mer nucleotides. There are \(4\) different options for each of these \(3\) positions {A, U, G, C}, leading to a total of \(64\) possible combinations. Additionally, 5 special tokens are added to the vocabulary: classifier token [CLS], separator token [SEP], unknown token [UNK], padding token [PAD], and masking token [MASK]. For finetuning, 64 codons, [CLS], [SEP] and [PAD] are used.

Since CodonBERT has an input length limit of \(1024\), we designed a hierarchical architecture in our CodonBERTER model (Figure ^[1]). Specifically, to accommodate the input length constraints of CodonBERT, the input mRNA sequence is first divided into \(n\) chunks, each with a maximum length of \(1022 \times 3\) nucleotides. These chunks are tokenized as codons, with [CLS] and [SEP] tokens appended to the start and end of each chunk, respectively. Unfilled chunks at the end of the sequence are padded with [PAD] tokens. After applying position and codon embeddings, chunks from the same mRNA sequence are stacked to form mini-batches, which are then fed into CodonBERT.

In CodonBERT, we use 12 layers of bidirectional transformer encoders, each with 12 self-attention heads, producing a hidden representation for each codon with a size of \(768\). The multi-head self-attention mechanism captures both long- and short-range interactions within the mRNA sequence, which influence translation efficiency. Each transformer encoder layer includes a feed-forward neural network that applies a non-linear transformation to the hidden representation obtained from the self-attention mechanism. Additionally, residual connections are applied around both the multi-head attention and feed-forward networks.

The output from CodonBERT is then flattened to form linearized codon embeddings for the input mRNA sequence. To bridge the CodonBERT output and the downstream model, we use an adaptor layer that applies a non-linear transformation to these embeddings. The transformed codon embeddings, combined with position information, are subsequently passed into a shallow transformer encoder block consisting of two transformer layers with eight self-attention heads. This block is designed to further extract both subtle long- and short-distance relationships in the sequence. Finally, protein expression levels are predicted using a dense layer that processes the extracted features.

Due to the lack of large-scale synonymous variant datasets of eukaryotes, we train, evaluate and test CodonBERTER (\(n=2\)) on the mRFP Expression dataset ^[6] profiling protein production levels for ~1500 gene variants in E. coli as a feasible preliminary attempt at this stage. We freeze the weights of attention layers in CodonBERT and fine-tune them using LoRA^[4], and train the other parameters of CodonBERTER from scratch. The training process takes 6h with a batch size of 12 around 100 epochs on a NVIDIA V100 PCIE 16GB.

In this section, we test several methods on the mRNA translation efficiency task, and finally compare them with our CodonBERTER model.

For a quick start, we select two typical predictors of mRNA expression, codon adaptation index (CAI) ^[8] and minimum free energy (MFE) as baseline, which are also used in LinearDesign algorithm ^[10]. CAI, defined as the geometric mean of the codon optimality of each codon in the mRNA \( \mathbf{r} \):

\[ \mathrm{CAI}(\mathbf{r}) = \sqrt[\frac{|\mathbf{r}|}{3}]{\prod_{0 \leq i < \frac{|\mathbf{r}|}{3}} w(\mathrm{codon}(\mathbf{r}, i))} \]

where \( \mathrm{codon}(\mathbf{r},i)=r_{3i}r_{3i+1}r_{3i+2} \) is the \(i\)th triplet codon in \( \mathbf{r} \), and \( w(c) \) is the relative adaptiveness of codon \( c \), defined as the frequency of \( c \) divided by the frequency of its most frequent synonymous codon \( 0< w(c) < 1 \).

\[ \mathrm{MFE}(\mathbf{r})=\min_{\mathbf{s}\in\mathrm{structures}(\mathbf{r})}\Delta G^\circ(\mathbf{r},\mathbf{s}) \]

where \( \mathrm{structures}(\mathbf{r}) \) is the set of all possible secondary structures for mRNA sequence \( r \), and \( \Delta G^\circ(\mathbf{r},\mathbf{s}) \) is the free-energy change of structure \( s \) for mRNA \( r \) according to an energy model.

The mRNA sequences in this dataset were divided into 3 groups according to CAI level, \(\mathrm{CAI_{high}}\) with overall medium CAI of 0.67, \(\mathrm{CAI_{medium}}\) 0.67 and \(\mathrm{CAI_{low}}\) 0.41. After normalizing the metrics to \([0,1]\), we calculate the Pearson and Spearman correlation coefficients between CAI/MFE and expression. Figure 2 shows the result. The first line shows the relationship between CAI/MFE of all mRNAs and protein expression levels. This is generally consistent with expectations, that is, higher CAI may mean faster translation, higher MFE may mean more longer half-life, and both are statistically significant correlated with higher protein expression. However, CAI and protein expression show only moderate positive correlation, while MFE show weak positive correlation, and the correlation with protein expression is worse, or even disappear, within the three CAI-level groups. Interestingly, the second graph in the first row seems to show that the MFEs of different CAI groups is not evenly distributed, suggesting that changes in CAI often lead to changes in MFE, and the two are not independent variables.

Next, we want to check if considering both CAI and MFE will improve the model's performance compared to using each variable independently. A multiple linear regression model (Figure 3, Table 1) taking both CAI and MFE as predictor variables was built. The Variance Inflation Factor (VIF) shows that there is no significant multicollinearity between the normalized CAI and MFE. However, the model has an adjusted \( R^2 \) of \( 0.126\), which means even when two variables are introduced, the explanatory power of the model does not increase and is even worse than that of the CAI single variable model. The negative correlation between normalized MFE and protein expression suggests that MFE may not bring more valuable information to the model but introduce irrelevant noise, which reduces the overall explanatory power of the multiple regression model.

\[ \begin{matrix} \begin{array}{llllllccc} \hline \text{ R-squared: } & 0.127& & \text{VIF}& \text{coef}&\text{std err}&\text{t}& P>|t|&[0.025&0.975] \\ \text{ Adj. R-squared: } & 0.126 & \text{ const } & 14.951396 & 0.6446 & 0.017 & 38.073 & 0 & 0.611 & 0.678 \\ \text{ F-statistic: } & 106.2 & \text{ CAI } & 1.410335 & 0.2424 & 0.018 & 13.677 & 0 & 0.208 & 0.277 \\ \text{ Prob (F-statistic): } & 8.99 \mathrm{E-44} & \text{ MFE } & 1.410335 & -0.1106 & 0.035 & -3.141 & 0.002 & -0.18 & -0.042 \\ \hline \end{array}\\ \text{Table 1: Summary of multiple linear regression model with CAI and MFE} \end{matrix} \]

Further, we want to see how the machine learning model performs in this task. In order to better learn sequence information, base-level features are embedded as input, including base pairing probability, one-hot coding of base types, and combination of the two. Two machine learning methods include Random Forest Regressor (RFR) and LASSO (Least Absolute Shrinkage and Selection Operator) are tested. As shown in the original paper, after embedding base level information, the accuracy of the model is greatly improved compared with the linear regression model based on RNA-level metrics. But the correlation coefficients between the prediction and the ground-truth of these methods are all \(< 0.8 \) (Figure 4) ^[6].

We ultimately employed CodonBERTER to address this task. After the model reached convergence (Figure 5A), we achieved a Spearman \( \rho \) of 0.87 and a Pearson \( r \) of 0.86 on the test dataset (Figure 5B), significantly outperforming other competing methods as reported in Table 2 ^[5]. The hierarchical model we designed for processing long mRNA sequences exhibits minimal performance loss compared to the original CodonBERT, suggesting that the downstream attention-based shallow transformer encoder block is highly effective in extracting and transferring sequence information. As shown in Figure 5C, the model's predictions are not strongly or significantly correlated with classical indicators like the CAI or MFE, suggesting that the model may have identified additional factors relevant to mRNA translation efficiency. Even when grouping the test data based on CAI levels, the model continues to exhibit strong predictive performance, although there is a slight reduction in accuracy for the high CAI group, which may be attributed to the limited number of training samples available for this category.

We also explored the interpretability of the model by analyzing attention weights from the final layer of the transformer decoder. After inputting an mRNA sequence, we observed that the model predominantly focuses on codons in the 10–20 range at the start of the sequence ^[3]. This region is typically associated with translation initiation, a rate-limiting step in protein synthesis, as well as with 5' end mRNA degradation, which affects mRNA stability. Therefore, the 5' end of the mRNA may play a key role in determining overall translation efficiency. Additionally, certain high attention weights located mid-sequence may reflect codon preferences or secondary structures that facilitate ribosome scanning, while peaks at the sequence's end could be associated with 3' end degradation processes.

\[ \begin{matrix} \begin{array}{ccc} & \text{Model} & \text{Spearman } \rho \\ \text{mRNA-based} & \text{CAI for single variable linear regression} & 0.44 \\ \text{Nucleotide-based} & \text{plain TextCNN} & 0.62 \\ & \text{RNABERT+TextCNN} & 0.40 \\ & \text{RFR} & 0.74 \\ & \text{LASSO} & 0.75 \\ & \text{RNA-FM+TextCNN} & 0.80 \\ \text{Codon-based} & \text{TF-IDF} & 0.57 \\ & \text{plain TextCNN} & 0.78 \\ & \text{Codon2vec+TextCNN} & 0.77 \\ & \text{CodonBERT} & \textbf{0.88} \\ & \text{CodonBERTER} & \textbf{0.87} \\ \end{array} \\ \text{Table 2: Model performance on mRFP translation efficiency prediction}^5 \end{matrix} \]

Building on the strong performance of CodonBERTER, we have developed a next-generation mRNA sequence design algorithm (Algorithm 1) that leverages the Codon Adaptation Index (CAI) for efficient search of solution spaces, combined with deep learning models for precise decision-making. In this approach, users input the desired protein amino acid sequence and target protein expression level. The algorithm then follows these steps:

1. (a) Initial Search Space Generation: A preliminary search space is generated based on a preset CAI value, which provides an estimate of codon usage that might influence translation efficiency.
2. (b) Prediction with CodonBERTER: The mRNA sequences from the initial search space are input into CodonBERTER to predict the corresponding protein expression levels. This step ensures that the model captures the complex relationships between mRNA sequences and translation efficiency, beyond what simple metrics like CAI can provide.
3. (c) Iterative Adjustment of CAI: If the predicted protein expression level does not meet the user's specified requirements, the CAI of the search space is adjusted accordingly: If the predicted expression is lower than the target, the CAI of the search space is increased to favor codons associated with faster translation. Conversely, if the predicted expression is higher than desired, the CAI is decreased to prioritize slower-translating codons.
4. (d) Iteration: The process is repeated, iteratively refining the search space and re-predicting expression levels using CodonBERTER, until an optimal mRNA sequence that satisfies the target expression level is identified.

This hybrid approach efficiently narrows down the solution space through CAI while leveraging CodonBERTER's deep learning capabilities for precise, context-aware predictions of protein expression, optimizing both speed and accuracy in mRNA sequence design.

In summary, we have modified an mRNA LLM to enable end-to-end prediction from mRNA sequence to protein expression level. Building upon this, we propose a next-generation mRNA sequence optimization algorithm. This work highlights the potential of deciphering the RNA language, which could revolutionize synthetic biology. However, it is important to recognize that the current RNA LLMs may face performance limitations when relying solely on unsupervised training. The model we trained on a dataset of single genes with uniform length is likely to suffer from inadequate generalization. Thus, there is an urgent need to develop high-throughput experimental techniques to generate high-quality datasets, which will be essential for fully exploring the potential of RNA LLMs.

1. Chen Bao, Mingyi Zhu, Inna Nykonchuk, Hironao Wakabayashi, David H. Mathews, and Dmitri N. Ermolenko. Specific length and structure rather than high thermodynamic stability enable regulatory mRNA stem-loops to pause translation. Nature Communications, 13, 2021.
2. Jiayang Chen, Zhihang Hu, Siqi Sun, Qingxiong Tan, Yixuan Wang, Qinze Yu, Licheng Zong, Liang Hong, Jin Xiao, Irwin King, et al. Interpretable RNA foundation model from unannotated data for highly accurate RNA structure and function predictions. arXiv preprint arXiv:2204.00300, 2022.
3. Wanjun Gu, Tong Zhou, and Claus O. Wilke. A universal trend of reduced mRNA stability near the translation-initiation site in prokaryotes and eukaryotes. PLoS Computational Biology, 6, 2010.
4. J. Edward Hu, Yelong Shen, Phillip Wallis, Zeyuan Allen-Zhu, Yuanzhi Li, Shean Wang, and Weizhu Chen. LORA: Low-rank adaptation of large language models. arXiv, abs/2106.09685, 2021.
5. Sizhen Li, Saeed Moayedpour, Ruijiang Li, Michael Bailey, Saleh Riahi, Lorenzo Kogler-Anele, Milad Miladi, Jacob Miner, Fabien Pertuy, Dinghai Zheng, Jun Wang, Akshay Balsubramani, Khang Tran, Minnie Zacharia, Monica Wu, Xiaobo Gu, Ryan Clinton, Carla Asquith, Joseph Skaleski, Lianne Boeglin, Sudha Chivukula, Anusha Dias, Tod Strugnell, Fernando Ulloa Montoya, Vikram Agarwal, Ziv Bar-Joseph, and Sven Jager. CodonBERT large language model for mRNA vaccines. Genome Research, 2024.
6. Thijs Nieuwkoop, Barbara R. Terlouw, Katherine G Stevens, Richard A. Scheltema, Dick de Ridder, John van der Oost, and Nico J. Claassens. Revealing determinants of translation efficiency via whole-gene codon randomization and machine learning. Nucleic Acids Research, 51:2363 – 2376, 2023.
7. Anastasiia Iu. Paremskaia, Anna A Kogan, Anastasiia Murashkina, Daria A. Naumova, Anakha Satish, I. S. Abramov, Sofya G. Feoktistova, Olga N Mityaeva, Andrei A. Deviatkin, and Pavel Yu Volchkov. Codon-optimization in gene therapy: promises, prospects and challenges. Frontiers in Bioengineering and Biotechnology, 12, 2024.
8. Paul M. Sharp and W. H. Li. The codon adaptation index – a measure of directional synonymous codon usage bias, and its potential applications. Nucleic Acids Research, 15(3):1281–95, 1987.
9. Qiushuang Wu, Santiago Gerardo Medina, Gopal Kushawah, Michelle Lynn DeVore, Luciana Andrea Castellano, Jacqelyn M. Hand, Matthew Wright, and Ariel Alejandro Bazzini. Translation affects mRNA stability in a codon-dependent manner in human cells. eLife, 8, 2019.
10. He Zhang, Liang Zhang, Ang Lin, Congcong Xu, Ziyu Li, Kaibo Liu, Boxiang Liu, Xiaopin Ma, Fanfan Zhao, Huiling Jiang, Chunxiu Chen, Haifa Shen, Hangwen Li, David H. Mathews, Yujian Zhang, and Liang Huang. Algorithm for optimized mRNA design improves stability and immunogenicity. Nature, 621:396 – 403, 2020.

The goal of using microneedles (MNs) as a delivery system for AAV vectors is to achieve efficient, targeted, and minimally invasive transfection of tissues. Upon insertion, the microneedles create microchannels in the epidermis, allowing for the efficient release of AAV into the deeper skin layers. The AAV is then free to diffuse through the interstitial fluid of the epidermis and into the target cells, where it can express therapeutic genes. Understanding the diffusion and transport dynamics of AAV particles after MN injection is critical for optimizing therapeutic efficacy¹.

AAV particles, with a size of around 25–30 nm, are small enough to diffuse effectively through tissues. The unique design of microneedles allows for localized delivery of these viral particles into the basal layer of the epidermis. This strategy reduces systemic side effects and enhances localized gene delivery.

We assume that the tissue is isotropic and homogeneous, meaning diffusion properties are uniform in all directions. The tissue thicknesses is of 15 μm for the stratum corneum and 100 μm for the viable epidermis ⁴. AAV nanoparticles are considered as 30nm nanoparticles in the diffusion model^2,3. The injection site is modeled as the source of AAV distribution, and diffusion is measured as time progresses. The evenly spaced microneedles on the support patch allow for the selection of a representative elementary volume (REV), as shown in Fig.2.

The core of the mathematical model revolves around Fick's laws of diffusion. We solve the governing equations for microneedle-mediated transdermal delivery of AAV using a 2D axisymmetric model.

AAV are assumed to be uniform before entering the skin tissue. Therefore, the AAV concentration in MN is set.

The concentration of AAV in the viable epidermis is driven by convection with interstitial fluid flow and diffusion due to concentration gradients. The governing equation is:

\[ \frac{\partial C_{AAV, VE}}{\partial t} = D_{AAV, VE} \nabla^2 C_{AAV, VE} \]

• \( C_{AAV, VE} \) is the AAV concentrationin the viable epidermis.
• \( D_{AAV, VE} \) is the AAV diffusivity in the viable epidermis.

AAV transport in the papillary dermis includes diffusion, convection, exchange with blood, and release:

\[ \frac{\partial C_{AAV, PD}}{\partial t} = D_{AAV, PD} \nabla^2 C_{AAV, PD} - \nabla \cdot (v_{is} C_{AAV, PD}) - Ex(C_{AAV, BL}, C_{AAV, PD}) \]

• \( C_{AAV, PD} \) is the AAV concentration in the papillary dermis.
• \( D_{AAV, PD} \) is the AAV diffusivity.
• \( Ex(C_{AAV, BL}, C_{AAV, PD}) \) is the exchange rate with blood.

After entering the blood, AAV can also continuously release the payload and be cleared by the organs such as the kidney. The concentration can be calculated by

\[ \frac{d C_{\mathrm{AAV}, \mathrm{BL}}}{d t}=\frac{V_{\mathrm{PD}} N}{V_{\mathrm{di}, \mathrm{AAV}}} E x\left(C_{\mathrm{AAV}, \mathrm{BL}}, C_{\mathrm{AAV}, \mathrm{PD}}\right)- k_{\mathrm{clr}, \mathrm{AAV}} C_{\mathrm{AAV}, \mathrm{BL}} \]

• \(k_{\mathrm{clr}, \mathrm{AAV}}\) is the plasma clearance rate.
• \( V_{\mathrm{PD}} \) is the volume of the papillary dermis surrounding each microneedle.
• \( V_{\mathrm{di}, \mathrm{AAV}}\) is the distribution volume of AAV.
• \(N \) is the number of microneedles in a patch.

Here are two parameter tables¹:

\[ \begin{array}{l} \text { Table 1. Model parameter for tissue properties. }\\ \begin{array}{lllll} \hline \text { Symbol } & \text { Parameter } & \text { Unit } & \text { VE } & \text { PD } \\ \hline \kappa_{\text {tis }} & \begin{array}{l}\text { Permeability to } \\\text { interstitial fluid }\end{array} & \mathrm{m}^{2} & 1.0\times 10^{-16} & 1.0\times 10^{-16} \\ \\ \mu_{\mathrm{is}} & \begin{array}{l} \text { Viscosity of interstitial } \\\text { fluid }\end{array} & \mathrm{Pa} \cdot \mathrm{~s} & 7.8\times 10^{-4} & 7.8 \times 10^{-4}\\ \\\rho_{\mathrm{is}} & \begin{array}{l}\text { Density of interstitial } \\\text { fluid }\end{array} & \mathrm{kg} / \mathrm{m}^{3} & 1000 & 1000 \\ \\\pi_{\mathrm{bl}} & \begin{array}{l}\text { Osmotic pressure of } \\\text { blood }\end{array} & \mathrm{Pa} & - & 2670 \\ \\\pi_{\mathrm{is}} & \begin{array}{l} \text { Osmotic pressure of } \\ \text { interstitial fluid } \end{array} & \mathrm{Pa} & - & 1330 \\ \\ \sigma_{\mathrm{T}} & \begin{array}{l} \text { Osmotic reflection } \\ \text { coefficient for blood } \end{array} & - & - & 0.91 \\ & \begin{array}{l} \text { proteins } \end{array} & & \\ \\ L_{\mathrm{bl}} & \begin{array}{l} \text { Hydraulic conductivity } \\ \text { of the blood vessel wall } \end{array} & \mathrm{m} / \mathrm{Pa} / \mathrm{s} & - & 2.7 \times 10^{-12}\\ \\ p_{\mathrm{bl}} & \text { Intracapillary pressure} & \mathrm{Pa} &-&2080 \\ \\ {S_{\mathrm{bl}}/ V_{\mathrm{tis}}}&\begin{array}{l} \text { Capillary surface area} \\ \text { per tissue volume } \end{array} &\mathrm{m^{-1}}&-&6.0\\ \hline \end{array} \end{array} \]

\[ \begin{array}{l} \text { Table 2. Model parameters for therapeutic agent properties. }\\ \begin{array}{llll} \hline \text { Symbol } & \text { Parameter } & \text { Unit } & \text { Nanoparticle } \\ \hline \\ D_{\mathrm{tis}} & \begin{array}{l} \text { Diffusivity in tissue of VE, PD } \\ \end{array} & \mathrm{m^2/s} & 1.0 \times 10^{-13}\\ \\ K & \begin{array}{l}\text {Partition coefficient between } \\\text { MN and tissue of VE, PD}\end{array} & - &1.0 \\ \\\sigma & \begin{array}{l}\text { Osmotic reflection coefficient} \\\text { for blood proteins}\end{array} & - & 1.0 \\ \\ P& \begin{array}{l} \text { Vascular permeability } \\ \end{array} & \mathrm{m/s}& 1.0 \times 10^{-9} \\ \\ C_{\mathrm{in}} & \begin{array}{l} \text { Administration dose} \\ \end{array}& \mathrm{M} & 1.0 \\ \\ k_{\mathrm{clr}} & \begin{array}{l} \text {Clearance rate in blood } \\ \end{array} & \mathrm{s^{-1}}& 5.0 \times 10^{-5}\\ \\ V_{\mathrm{dis}} & \begin{array}{l} \text {Distribution volume} \end{array} & \mathrm{m^3} & 1.8 \times 10^{-2}\\ \hline \end{array} \end{array} \]

Based on the simulation results, the use of 115 nm long microneedles for AAV vector injection leads to a low concentration in the papillary dermis (PD), with the majority of the AAV particles remaining concentrated in the upper viable epidermis (VE). This finding suggests that most AAV vectors remained in the upper layer of the epidermis, potentially due to the rapid turnover of this skin layer. The constant regeneration of the epidermis may lead to a shorter effective period for a single injection, as newly formed cells could quickly replace those containing the AAV. Consequently, the therapeutic effects may diminish rapidly, necessitating multiple injections or alternative delivery methods to maintain sustained levels of AAV within the target tissue.

Due to the lack of in vivo experimental data for AAV, such as AAV transduction efficiency and XPC expression levels, it is currently not possible to determine how to design the microneedles or what dose to inject. Further experimental research is still needed.

Ultimately, understanding these dynamics could inform future designs for microneedle systems, allowing for more effective delivery of AAV vectors, maximizing therapeutic potential while minimizing systemic side effects.

1. Newell B, Zhan W. Mathematical modelling of microneedle-mediated transdermal delivery of drug nanocarriers into skin tissue and circulatory system. J Controlled Release. 2023;360:447-467. doi:10.1016/j.jconrel.2023.07.011
2. Microvesicle-associated AAV Vector as a Novel Gene Delivery System. Mol Ther. 2012;20(5):960-971. doi:10.1038/mt.2011.303
3. Kato M, Ishikawa S, Shen Q, et al. In situ-formable, dynamic crosslinked poly(ethylene glycol) carrier for localized adeno-associated virus infection and reduced off-target effects. Commun Biol. 2023;6(1). doi:10.1038/s42003-023-04851-w
4. Predicting viable skin concentration: Diffusional and convective drug transport. J Pharm Sci. 2021;110(7):2823-2832. doi:10.1016/j.xphs.2021.03.012