We consider single 20μm silicon photodiodes featuring an inner electrode (with either a square or fractal geometry) and an outer, grounded electrode (Figure 1). Both electrodes are 250 nm tall and are composed of titanium nitride (TiN), a commonly used retinal implant electrode material (Zrenner et al., 2011; Stingl et al., 2015). The silicon area (Figure 1, green) is 16 × 16μm and is surrounded by a 500 nm wide insulating layer (Figure 1, yellow). The bounding area (Figure 1, dashed white lines) for the square electrodes is varied between 50 and 200 μm². The construction of the fractal electrodes is as follows.

Mathematically exact fractals can be constructed by scaling an initial seed pattern and then iterating the scaled pattern toward increasingly fine size scales. The scaling rate, L, is set by the number of new patterns created, N, and D, according to the equation

where 1 ≤ D ≤ 2. Throughout this paper we model branched “H-tree” fractal electrodes. Figure 2 illustrates H-tree fractals which hold D fixed at 2.0 and increases the iterations from 1 to 2 to 3, and also H-trees which hold the iterations fixed at 3 and increases D from 1.4 to 1.7. In general, the H-tree electrode becomes more space filling for increasing iterations and increasing D. Each fractal electrode features line widths of 160 nm and a fixed bounding area of 15.4 × 15.2 μm. This line width was selected due to its ease of fabrication and also because it prevents the branches at different iterations from overlapping. In total, 13 electrode geometries were studied: 4 square electrodes with covering areas of 50, 100, 150, and 200 μm² and 9 fractal electrodes from each combination of D values of 1.4, 1.7, and 2.0 and iterations of 1, 2, and 3. The covering area (i.e., the area of the top surface of each electrode) and total surface area (i.e., including the sidewall surfaces) of each electrode is given in Table 2.

The general strategy applied to the 3 simulation steps used to determine electrode, neuron, and photodiode responses is to mesh 3-dimensional geometries into a set of nodes, establish an equivalent circuit model between nodes (e.g., Figure 3 for 2-dimensional illustrations), and calculate the node voltages using modified nodal analysis (MNA) (Ho et al., 1975). Briefly, MNA determines node voltages by applying Kirchhoff's current conservation rule at each node along with the appropriate boundary conditions. For n node voltages, V→  = (V₁, …, V_n), and m applied voltage sources, V→ ^app =  (V1app,…,Vmapp), the MNA system of equations is given by

where G is an n × n matrix containing conductance elements between nodes, A is an m × n matrix that sets boundary conditions to the applied voltages and only contains zeros and ones, I→ = (I₁, …, I_m), gives the m currents flowing through the applied voltage sources, and I→ ^app =  (I1app,…,Imapp) applies current sources to the n nodes. The lower right m × m matrix is zero. The system of equations is solved using the package SuperLU (Demmel, 1999; Li and Demmel, 2003).

The MNA algorithm outlined above is used to characterize the current and voltage generated by each photodiode under illumination. The photodiode is first recreated as a 2-layer cubic mesh featuring TiN electrode nodes in the top layer and semiconducting silicon nodes in the bottom layer (Figure 3A). The node-to-node impedances feature an electrode resistance between metal nodes, a sheet resistance between semiconducting nodes, and a contact resistance between metal and semiconducting nodes. The equivalent circuit model also includes the load impedance magnitude, |Z|, between the inner and outer electrodes (see section 2.4 and Equation 5). Under illumination, the photodiode current is modeled as an array of current sources (i.e., photocurrents generated from the incident radiation) in parallel with diodes (i.e., “dark” currents—thermal recombination between the electron and hole charge carriers that occurs even in the absence of light) (Figure 3A). Photocurrents are included only for nodes which are exposed to the radiation (i.e., not blocked by the inner electrode). The net current for node j, Ijapp, is given by an ideal diode under illumination according to

where a is the node's top surface area, I_sc is the short-circuit current (current when |Z| = 0), I_dark the “dark” reverse current, I_rad the irradiance (radiant power per unit area), R the photodiode responsivity (amps generated per incoming watt of radiation), and J₀ the dark current density at 0 V. The dark current density is estimated to be 100 nA/cm² by comparing to similar microphotodiode subretinal implants (Chow et al., 2004; Wang et al., 2012). The thermal voltage, V_T = 0.0268 V, is the value at the body's temperature of 310 K. V_j is the voltage at node j. Semiconducting nodes below the top-contact only feature a dark current. The only quantity inserted in V→ ^app (Equation 2) is setting the ground potential to 0 V. For each photodiode, the MNA equation is solved iteratively using a global Newton method (Bank and Rose, 1981; Ceric, 2005) to determine the node voltages, V→, and the current flowing through the load impedance, I. We note that each electrode exhibited a maximum variation in metal node voltages of less than 1.3%. Each inner electrode is therefore approximately equipotential, and is measured by a single voltage at the central metal node, termed V. The relationship between I and V, gives the photodiode's IV curve.

In addition to the IV curves, one common characterization of subretinal implant photodiodes that occurs on the laboratory benchtop, prior to implantation, is to measure the open-circuit voltage, V_oc, by leaving the connection between the inner and outer electrodes open (i.e., infinite |Z|). For a given irradiance, V_oc only depends on the photodiode parameters and the inner electrode geometry. Here, the open circuit voltage, V_oc, can be estimated by

where A_tot is the total photodiode area and A_pd is the photodiode area not blocked by the electrode.

We previously described how the MNA algorithm outlined in section 2.2 is used for calculating the extracellular voltages generated by the electrodes (see Watterson et al., 2017 for a detailed explanation). Briefly, a 1 mm³ cubic domain containing the inner electrode, the outer grounded electrode, and the extracellular space is meshed into a set of tetrahedral nodes. Next, an equivalent circuit model is created which defines the node-to-node impedances. The fluid-fluid nodes are resistive (R_f), while the fluid-electrode nodes feature a capacitor (C_dl) and resistor (R_ct) in parallel, which model charge screening and reversible oxidation-reduction reactions at the electrode, respectively (Figure 3B) (Merrill et al., 2005). The fluid resistivity is taken to be 3,500 Ω cm (Kasi et al., 2011). The applied voltage boundary conditions, V→ ^app in Equation (2) are set as follows. An oscillating voltage, V_e =Ve^2πift, (where the value of V is inputted from the photodiode simulations) is applied to the inner electrode while the outer electrode is held at 0 V. The remaining boundary conditions are set to be insulating for the plane in which the electrode is located, and 0 V at the other 5 faces of the cubic domain. There are no applied current sources in this portion of the simulations so I→ ^app =  0→ (Equation 2). Having established the equivalent circuit model along with the boundary conditions, Equation (2) can be solved for the n complex valued node voltages and m complex currents through the boundary condition nodes. The voltage in the electrolytic fluid outside the neuron is termed V_f. The load impedance magnitude, |Z|, (which is set by the network of R_f, C_dl, and R_ct components) can also be calculated by

where |I| is the current leaving the inner electrode. Additionally, we calculate the charge density, Q_ph, at each node delivered on the electrode surface per positive phase of voltage by

The neuron stimulation is described in detail elsewhere (Watterson et al., 2017). Briefly, the extracellular voltages, V_f, calculated in section 2.4 induce a change in the membrane potential, ΔV_m = ΔV_i − ΔV_f, in the bipolar neurons located near the electrode, where ΔV_i is the change in internal potential of the neuron and ΔV_f is the change in electrolytic fluid potential. In turn, these bipolar neurons pass their signal downstream to retinal ganglion cells when ΔV_m reaches a minimum of 15 mV at the bipolar neuron's soma (Yang and Wu, 1997).

Our model bipolar neurons are 100 μm long with a 10 μm soma centered 30 μm above the electrode surface, and dendrites which are just above the electrode's top surface (Wassle et al., 2009; Masland, 2012). In our simulations, each neuron features a cubic mesh. Passive rod bipolar neurons are quantified by a membrane capacitance of 1.1 μF/cm² in parallel with a membrane resistance of 2.4 × 104 Ωcm², along with an internal cytoplasmic resistivity of 130 Ω cm (Oltedal et al., 2009). For the applied stimulation frequencies used here (1 kHz), the resistive impedance is more than 2 orders of magnitude higher than the capacitive impedance. We therefore ignore the resistive component and create an equivalent circuit model containing solely membrane capacitances and internal cytoplasmic resistances (Figure 3C). The real and imaginary parts of the extracellular voltages obtained in section 2.4 are mapped onto the outside of the neuron's membrane and serve as a set of applied voltage sources, V→ ^app in Equation (2). The MNA equation is then solved to obtain the neuron's internal voltage at each node.

We first consider the load impedance, |Z|, because this will determine how close the photodiode operates at to open or closed-circuit. The load impedance for the square electrode is found to decrease with increasing electrode size (Figure 4). This is expected because the geometric contribution to the load impedance is inversely proportional to the inner electrode's effective surface area and directly proportional to the distance between the inner and outer electrodes. The fractal electrode reduces its impedance relative to the square by increasing its effective surface area (by maximizing the surface area via the large number of branch sidewalls) and decreasing the separation between inner and outer electrodes. This leads to a general trend of decreasing impedance for increasing D value and increasing iterations (Figure 4).

Figures 5A–C show the IV curves for photodiodes with square and fractal electrode geometries under an illumination of I_rad = 10 mW/mm². Each electrode exhibits an open-circuit voltage of ~0.40 V and a short-circuit current proportional to the exposed photodiode area. In a conventional solar cell, the load impedance would be chosen to maximize the power generated. However, here each photodiode has a load impedance set by the electrode geometry. The black dot on each trace in Figures 5A–C shows the operating point on the IV curve set by the impedances reported in Figure 4.

The voltages generated by each electrode geometry as a function of irradiance display several common characteristics (Figures 5D–F). Firstly, at low voltages the slope of each trace is given by ΔV/ΔI_rad =R|Z|A_pd, where A_pd is the unblocked photodiode area. The electrodes with smaller covering areas have both large |Z| and A_pd and therefore generate relatively high voltages at the lower intensities. Secondly, as the voltage begins to approach the open circuit voltage, increasing illumination intensity provides minimal increases in the electrode voltage.

The results of section 3.1 highlight the importance of electrode geometry when determining the voltage generated for a given illumination. However, electrode geometry also influences how the field from this voltage extends into the extracellular liquid and this can lead to competing considerations. For instance, Figures 6A,D show the effect of increasing a square electrode's area from 100 to 150 μm². As expected from the decrease in |Z| and A_pd, the larger electrode's voltage decreases significantly and the field does not therefore extend as far vertically into the liquid as the smaller electrode's field. However, the field from the larger electrode has the advantage of extending further horizontally within the pixel. An inevitable consequence of the square design therefore is that fields that extend far vertically do not extend far horizontally and vice versa.

The fractal design offers a potential solution for optimizing this competition. The fractal electrode generates high voltages for a given illumination (Figure 5). Furthermore, its maximal capacitance (due to the large surface area generated by the branch sidewalls) allows a large amount of charge to reside on the electrode and this generates a large field for a given applied voltage, which will penetrate far vertically into the liquid. Because the electrode spreads further laterally than a square electrode for the same covering area, the fractal electrode's field will also extend far horizontally within the confined area of a single photodiode. We note that the horizontal spread does not significantly spread into the neighboring pixel due to the outer ground, and therefore the electrical crosstalk remains minimal (see Watterson et al., 2017). However, the presence of the gaps in the fractal design needs to be taken into account. Figures 6B,E show the fields for the D = 1.4 electrode; both the 2 and 3 iteration electrodes feature large gaps which reduce the extracellular voltage in the central region. For the D = 2.0 fractals shown in Figures 6C,F, increasing the number of iterations from 2 to 3 reduces the voltage but the field spreads out relatively uniformly across the entire pixel. Given that larger extracellular fields generally induce large depolarizations, ΔV_m, of the bipolar neurons, it is clear from the above that careful geometric optimization will be required to supply a large voltage which extends into the most extracellular space.

The stimulation efficiency for each design is determined by measuring ΔV_m, for a patch of 9 bipolar neurons directly above each electrode. Figure 7 depicts ΔV_m for a patch of 4 of the 9 neighboring bipolar neurons above electrodes under equivalent illuminations of 10 mW/mm². Because the 150 μm² square electrode (Figure 7A) blocks a larger percentage of the underlying photodiode and therefore has a lower voltage on the inner electrode, the neurons above the square depolarize less compared to the 2 iteration D = 1.4 and D = 2.0 H-trees (Figures 7B,C). Additionally, the fractal electrode's D value influences the field distribution in the extracellular space, leading to varying neural depolarizations. For instance, although the voltage on the 2 iteration D = 1.4 H-tree is slightly larger than the voltage on the 2 iteration D = 2.0 H-tree (0.38 vs. 0.36 V), the depolarizations are larger for neurons above the 2 iteration D = 2.0 H-tree (Figures 7B,C).

To quantify the stimulation efficiency, we define the electrode threshold stimulating voltage, V_thresh, as the electrode voltage at which all 9 neighboring bipolar neurons reach a somatic depolarization of ΔV_m = 15 mV. Previous experiments show this 15 mV condition results in stimulation of the downstream ganglion neurons Yang and Wu (1997). For square electrodes, increasing the electrode area reduces V_thresh due to an increase in capacitance. Likewise, increasing the capacitance for fractal electrodes either by increasing the number of iterations or increasing the D value leads a to lower V_thresh (Figure 8A). However, as discussed in section 3.1, increasing the electrode's covering area also reduces the voltage generated on the inner electrode. Therefore, efficient stimulation requires a careful optimization of supplying enough voltage from the photodiode and maintaining a low V_thresh.

Across all of the electrode patterns, the 2 iteration D = 2.0 fractal provides the best balance between these 2 competing factors (Figure 8). In particular, the incident radiation required to stimulate all neighboring bipolar neurons is 290% more for the best square electrode of 200 μm² than for the the 2 iteration D = 2.0 fractal. We note that at their threshold voltages, the maximum surface charge density, Q_ph, of the optimized electrodes are Q_ph = 0.67 mC/cm² for the 200 μ m² square electrode and Q_ph = 0.93 mC/cm² for the 2 iteration D = 2.0 H-tree. These charge densities are less than the 1 mC/cm² safety limit for TiN electrodes based on the charge densities that induce hydrolysis (Weiland et al., 2002).

So far, we have considered stimulating pulses operating at a frequency of 1 kHz. However, conventional subretinal implants being developed today use stimulating frequencies ranging from 250 Hz to 2 kHz (Zrenner et al., 2011; Mathieson et al., 2012; Lorach et al., 2015). In order to verify the fractal maintains a lower threshold irradiance at lower stimulating frequencies than 1 kHz, we repeated the above analysis for the 150 μm² square, the 200 μm² square, and the 2 iteration D = 2.0 fractal at a stimulating frequency of 250 Hz.

First, lowering the stimulating frequency causes a rise in the load impedance, |Z|, for each geometry due to an increase in capacitive impedance at the electrode-electrolyte interface. This increased |Z| leads to a larger voltage generated on the inner electrode (i.e., the operating point on the IV curve shifts to a higher voltage). For example, under 10 mW/mm² illumination, reducing the frequency from 1 kHz to 250 Hz causes an increase in the inner electrode voltage from 0.09 to 0.11 V for the 200 μm² square and from 0.34 to 0.36 V for the 2 iteration D = 2.0 fractal. Simultaneously though, the increased impedance leads to a smaller spreading in the extracellular field generated by each electrode (Figure 9). Additionally, the lower frequency causes smaller depolarizations in the bipolar neurons due to a higher capacitive membrane impedance.

Combining all of these factors, we find the threshold irradiances, I_thresh, necessary to depolarize all 9 surrounding neurons at 250 Hz are 90 mW/mm² for the 150 μm² square, 42 mW/mm² for the 200 μm² square, and 15 mW/mm² for the 2 iteration D = 2.0 fractal. Therefore, lowering the stimulating frequency from 1 to 250 kHz causes a reduction in I_thresh for the 200 μm² square and an increase in I_thresh for the 2 iteration D = 2.0 fractal. However, the fractal implant still requires 64% less irradiance intensity to stimulate all surrounding neurons than the best square design.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.